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ABSTRACT 


The  application  of  analysis  lattice  filters  to  the 
problem  of  determining  the  input  to  a  system  from 
observations  of  the  system’s  output  (i.e.,  deconvolution)  is 
discussed.  Both  linear  and  nonlinear  systems  are 
considered.  Lattice  filter  modeling  algorithms  (Levinson 
and  Schur)  are  presented. 

The  theory  of  least-squares  inverse  filters  is  reviewed. 
This  leads  to  a  discussion  of  the  lattice  filter,  which  in 
turn  leads  to  the  Generalized  Lattice  Theory.  The 
Generalized  Lattice  Theory  is  then  used  to  develop  a 
nonlinear  lattice  structure.  Simulations  show  that  the 
nonlinear  lattice  is  an  effective  inverse  filter  for  both 
linear  and  nonlinear  systems. 
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I .  INTRODUCTION 


I^he  problem  of  estimating  a  signal  based  upon 
observations  of  a  related  signal  is  one  of  the  most 
important  operations  in  signal  processing.  The  input  and 
output  signals  of  a  linear  system  are  related  by  the 
convolution  operation 

y  ( t )  =  h  ( t )  *  x ( t )  =f  h(t-T)x(T)d'f  (1.1) 

where  h(t)  is  a  causal,  linear  time-invariant  (LTI),  system 
impulse  response.  Since  this  thesis  deals  primarily  with 
discrete  digital  signals,  continuous  time  signals  are 
sampled  at  uniform  intervals,  T,  and  are  represented  by 
discrete  sequences  x ( nT )  =  x(n)  for  n  =  0,1,2,...,N. 

Discrete  convolution  for  a  LTI  system  is  defined  by 

n 

y(n)  =  h(n)  *  x(n)  =  £  h(n-m)x(m)  (1.2) 

m=-co 

and  is  shown  in  Figure  1.1.  The  basic  deconvolution 
problem  is  to  estimate  the  signal  x(n),  assuming  that  both 
y(n)  and  h(n)  are  known.  Figure  1.2  depicts  the  inverse 
filtering  process  of  recovering  the  input  signal  from  the 
output  signal  by  removing  the  system’s  impulse  response. 

Deconvolution  has  important  applications  in  a  variety 
of  fields:  Radar,  communications,  image  processing,  speech 


Observed 


Original 


Tf  'J 


synthesis,  seismology.  For  example,  in  image  processing, 
deconvolution  is  used  to  recover  the  representation  of  the 
original  object,  x(m,n),  from  the  representation  of  its 
image,  y(m,n),  by  removing  the  blurring  caused  by  the  opti¬ 
cal  system’s  point-spread  function,  h(m,n).  A  common 
problem  which  arises  in  geophysics  involves  the  deconvolu¬ 
tion  of  a  seismic  trace,  y(n),  into  the  approximately  known 
impulsive  waveform,  h(n),  and  the  desired  reflection 
response,  x(n),  which  reveals  the  structure  of  the  layered 
Earth.  In  other  applications,  h(n)  may  represent  the 
impulse  response  of  a  transmission  channel,  magnetic 
recording  medium,  or  measurement  device  which  broadens  and 
smears  (intersymbol  interference)  the  desired  message  x(n). 

There  have  been  numerous  approaches  to  the  linear 
deconvolution  problem,  including  least-squares  filtering, 
linear  inverse  theory,  linear  programming,  and  homomorphic 
signal  processing.  The  first  part  of  this  thesis  deals  with 
least-squares  filtering  techniques;  the  application  of 
Kalman,  waveshaping,  and  lattice  filters  to  deconvolution  is 
reviewed.  Next,  the  lattice  filter  discussion  is  extended 
to  the  theory  of  the  generalized  lattice  filter.  Finally, 
nonlinear  system  theory  is  briefly  reviewed,  and  the 
nonlinear  lattice  filter  is  developed  and  applied  to  the 
inverse  filtering  problem.  Computer  simulation  results  for 
both  the  linear  and  nonlinear  lattice  filters  are  presented. 
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A. 


INTRODUCTION 


This  chapter  begins  with  a  discussion  of  the  principles 
of  least-squares  filtering  theory  since  deconvolution  is 
essentially  an  inverse  filtering  process.  Three  specific 
types  of  least-squares  filters  are  then  introduced:  Kalman, 
waveshaping,  and  lattice  filters.  The  application  of  each 
of  these  three  filter  types  to  the  problem  of  linear  decon¬ 
volution  is  studied.  Finally,  the  chapter  concludes  with 
the  presentation  of  computer  simulation  results  obtained  for 
deconvolution  experiments  using  the  lattice  filter. 

The  goal  of  deconvolution  is  to  recover  the  input 
signal,  x(n),  to  a  system  based  upon  observed  values  of  the 
system’s  output  signal  yin).  An  optimal  processor  must  be 
determined  to  produce  the  best  possible  estimate  of  x(n) 
based  upon  present  and  past  values  of  the  output  y(n).  A 
traditional  measure  for  defining  the  "best"  signal  processor 
is  the  minimum  mean-squared  error  criterion.  Using  this 


criterion , 

the 

estimate  x(n) 

is  defined  by 

a  linear 

combination 

o  f 

the  observed 

values  y ( n  )  . 

Assuming 

causality,  and  that  the  observed  signal  is  windowed  to 
include  only  the  M  past  samples,  x(n)  is  written  as 


x(n)  =  Y  h(n,i)y(i)  ,  n  ^  0 
i  =  n-M 

where  h(n,i)  =  0  for  n  <  i.  The  estimation  error 
e(n)  =  x(n)  -  x(n).  To  find  the  optimum  signal 
the  h(n,i)  coefficients  which  minimize  the 
estimation  error  J,  where 


(2.1) 

is  given  by 
processor , 
mean-square 


2 

J  =  E [ e  (n) ] 


E  [  (  x  (  n  ) 


x  (  n  )  ) 


2 

]  , 


(2.2) 


must  be  determined.  This  is  known  as  the  Wiener  filter 
formulation  of  the  problem.  [Ref.  1 : pp .  113,116] 

The  least-squares  theory  of  filtering  began  in  the 
1940’s  with  the  work  of  Norbert  WIENER  [Ref.  2 : pp .  147-148]. 
WIENER  developed  a  frequency  domain  procedure  to  design 
optimum  filters,  where  optimality  was  defined  by  minimizing 
a  mean-square  error  performance  criterion.  The  Wiener 
filter  is  conventionally  applied  to  linear  time-invariant 
systems  with  stationary  statistics  when  it  is  desired  to 
separate  one  signal  from  another.  In  the  early  1950’s,  the 
Wiener  filter  was  extended  to  include  time  varying  and 
nonstationary  statistics,  but  the  calculations  are 
cumbersome  [Ref.  3:p.  1]. 

The  mean-square  estimation  error  J(n)  is  minimized  by 
setting  its  partial  derivatives,  with  respect  to  each  of  the 
filter  coefficients  h(n,i),  equal  to  zero. 


1  1 


Thus  : 


0h(n,  i  ) 


0 


(2.3) 


for  i  =  n-M,n-M+l . n.  This  yields  a  set  of  M+l  linear 

simultaneous  equations,  called  ths  orthogonality  equations 
where 


R  ( n , i )  =  E[ e ( n ) y ( i ) ]  =  0  (2.4) 

ey 

for  n-M  <_  i  <_  n  and  n  >_  0 .  R  (n,i)  is  the 

ey 

crosscorrelation  function  between  the  error  signal,  e(n), 
and  the  data  y(n).  Substituting  e(n)  =  x(n)  -  x(n)  into  the 
orthogonality  equations  gives  the  normal  ,  or  Wiener-Hopf 
equations : 

n 

E [ x ( n ) y ( i ) ]  =  Y  h ( n , k ) E [ y ( k ) y ( i ) ]  (2.5a) 

k=n-M 


or 

n 

R  (n,i)  =  £  h(n,k)R  (k,i)  .  (2.5b) 

xy  k=n-M  yy 


Since  the  autocorrelation  function  R  of  the  input  signal 

yy 

and  the  crosscorrelation  function  R  of  the  desired  output 

xy 

signal  with  the  input  signal  are  known  quantities,  the  M+l 
equations  can  be  solved  for  the  optimal  filter  weights 
h(n,i),  i=n-M,...,n.  If  data  vectors  are  defined  so  that 

T 

x(n)  =  [x(n-M),  x ( n-M+ 1 ),..., x ( n ) ]  (2.6a) 
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and 


y(n)  =  [y(n-M),  y(n-M+l) . y(n)]  (2.6b) 

then  the  M+l  Wiener-Hopf  equations  can  be  written  as 
T  T 

E[x(n)y(n)]=hE[y(n)y(n)]  (2.7) 

where  h  =  [h(n,n-M),  h ( n-M+ 1 ) , . . . , h ( n , n ) ] .  Assuming  that 
the  signals  x(n)  and  y(n)  are  stationary,  then  the  filter 
coefficients  are  time-invariant  and  h(n,k)  =  h(n-k);  and 

equation  (2.7)  can  be  written  in  matrix  form  as 


(2.8) 


Ryy(O) 

Ryy (  1 ) 

Ryy ( 2  )  ... 

Ryy ( M ) 

h  (  0 ) 

- 

Rxy ( 0 ) 

Ryy ( 1 ) 

Ryy ( 0 ) 

Ryy ( 1  )  .  .  . 

h  (  1  ) 

Rxy ( 1 ) 

Ryy ( 2 ) 

• 

Ryy (  1 ) 

Ryy ( 0  )  ... 

•  • 

h  (  2  ) 

• 

= 

Rxy ( 2 ) 

• 

• 

• 

Ryy ( M ) 

• 

• 

•  • 

• 

Ryy  ( 0 ) 

• 

• 

h(M) 

Rxy(M) 

Now,  the  optimal  coefficients  can  be  solved  by  inverting 
the  autocorrelation  matrix,  or  by  exploiting  the  matrix’s 
Toeplitz  structure  (all  elements  are  the  same  on  any  given 
diagonal)  to  employ  the  more  efficient  Levinson  algorithm. 
(Levinson’s  algorithm  will  be  discussed  in  the  development 
of  the  analysis  lattice  filter.) 

The  minimized  value  of  the  mean-square  estimation  error 
can  now  be  computed,  and  is  found  to  be 


2  n 

J  (n)  =  E[x  ( n ) 1  -  £  h ( n , i ) E [ y ( i ) x ( n ) ] 

min  i=n-M 


(2.9) 


>■ -V  •“•"t*  ■  .■'\0vv 


w 


:  W 


% 


=  R  (  0  )  -E[  x(  n  )  x  (n>]E(y(n>y  C  n  )  )  E[y(n)\<n) 


[Ref.  4:p.  1481.  These  values  correspoi.  i  to  the  diagonal 

elements  of  the  covariance  matrix  of  the  estimation  error, 


R  =  E  [  e  (  n  )  e  (n)] 


(2.10) 


where  e(n)  =  x(n)  -  x(n)  [Ref.  l:p.  120).  Furthermore,  the 


estimate  of  x(n)  is  given  by 


T  -1 


x  (  n  )  =  E  [  x  (  n  )  v  (  n  )  ]  E  [  £  (  n  )  v  (  n  )  1  y  (  n  ) 


=  hv ( n ) 


(2.11) 


This  estimate  can  be  thought  of  as  the  pro  'ection  of  the 
desired  signal  x(n)  onto  the  space  spanned  by  the  components 
of  the  observation  vector  y(n).  The  minimized  estimation 
error  vector  is  orthogonal,  or  normal,  to  the  estimate 
x(n).  [Ref.  4:p.  147] 

This  completes  the  overview  of  least-squares  filtering. 
The  following  sections  present  several  linear  deconvolution 
techniques  which  employ  this  criterion:  Kalman  filtering, 
spiking  filters,  and  lattice  filters. 


B.  KALMAN  FILTER 

In  the  early  1960’s,  R.E.  KALMAN  introduced  an  optimal 
recursive  filter  based  on  state-space  time-domain  methods 
[Ref.  2:pp.  267-268].  The  Kalman  filter  estimates  the  state 


of  a  linear  system,  and  is  optimal  in  the  sense  that  it  min¬ 
imizes  the  mean-square  error  of  the  state  estimate.  The 
Kalman  filter  is  useful  when  the  system  is  defined  by  state 
space  equations:  The  system  signals  are  represented  by 
random  processes,  and  the  data  observations  are 
contaminated  by  noise.  The  Kalman  filter  algorithm  processes 
measurement  data,  and  requires  a  priori  state  space  models 
(known  or  assumed)  of  the  system  and  measurement  dynamics. 
Also,  the  statistics  of  the  system  input  and  measurement 
noises,  as  well  as  initial  condition  information,  are 
required  to  produce  the  state  estimate.  This  process  is 
depicted  in  Figure  2.1.  Here,  the  discrete  Kalman  filter 
equations  will  be  presented,  and  then  their  application  to 
deconvolution  will  be  described. 

The  state  space  representation  of  the  discrete  system 
and  measurement  models  (see  Figure  2.2)  are  written  as  [Ref. 
2 : pp.  195-200]: 

x ( k )  =  F ( k- 1 ) x ( k- 1 )  +  w(k-l)  (2.12) 

z(k)  =  H(k)x(k)  +  v(k)  (2.13) 

where 


x  ( k ) 

=  ( n 

X 

1  ) 

system  state 

vector 

F  ( k ) 

=  (n 

X 

n ) 

transition  matrix 

w(k) 

=  ( n 

X 

1  ) 

system  noise 

vector 

z  ( k  ) 

=  ( m 

X 

1  ) 

measurement 

vector 

H(k)  =  (m  x  n)  observation  matrix 

v(k)  =  (m  x  1)  measurement  noise  vector. 

(Note  that  the  time  index  k  is  used  here  to  be  consistent 

with  the  discrete  Kalman  filter  literature.) 

The  noise  vectors  w(k)  and  v(k)  are  assumed  to  have  zero 

mean,  white,  Gaussian  distributions  with  covariance  matrices 

T  T 

of  Q(k)  =  E[w(k)w  (k)]  and  R(k)  =  E[v(k)v  ( k ) ] , 

respectively.  Additionally,  w  and  v  are  uncorrelated  so 

T 

that  E[w(k)v  ( j ) ]  =  0  for  all  k  and  j .  The  state  estimation 
error  is  defined  by 

e(ktk-l)  =  x(k)  -  x(k|k-l)  (2.14) 

and  the  associated  (n  x  n)  error  covariance  matrix  is 

T 

P<k|k-1)  =  E [ e ( k | k- 1 ) e  ( k | k-  1  >  ]  .  (2.15) 

An  updated,  a  posteriori  estimate  of  x(k)  is  obtained  from 
the  measurement  z(k)  and  the  a  priori  state  estimate 
x ( k | k- 1 )  by 

x ( k )  =  x(k|k-l)  +  K(k)[z(k)-H(k)x(k|k-1)  )  (2.16) 

where  K(k)  is  the  (n  x  m)  Kalman  gain  matrix.  The  a 
posteriori  error  is  given  by 

e(k)=x(k)-x(k)  (2.17) 


and  the  associated  a  posteriori  error 


covariance  matrix  is 


T 

P(k)  =  E  [  e  ( k  )  e  (k) ] .  (2.18) 

The  mean-square  estimation  error  criterion  is  minimized 
when 

T  T  -1 

K( k ) =P( k | k-1 )H  (k) [H(k)P(k|k-l)H  (k)+R(k)]  (2.19) 

This  optimal  value  of  the  Kalman  gain  matrix  minimizes  the 
individual  terms  along  the  main  diagonal  of  P(k). 
Substituting  the  optimal  gain  matrix  K(k)  into  the 
expression  for  P(k)  results  in 

P(k)  =  (I  -  K ( k ) H ( k ) )  P<k|k-1)  ,  (2.20) 

where  I_  is  the  identity  matrix.  In  order  to  compute 
equations  (2.16)  and  (2.19)  recursively,  the  a  priori 
estimates  x(k+l|k)  and  P(k+l|k)  must  be  determined  at  time 
k.  The  a  priori  estimates  are  given  by 

x ( k+ 1 | k )  =  F ( k ) x ( k )  (2.21) 

T 

P ( k+ 1 | k )  =  E [ e ( k+ 1 | k ) e  (k+llk)]  (2.22) 

T 

=  E[F(k)e(k)+w(k) ) ( F ( k ) e ( k ) +w( k ) )  ] 

T 

=  F(k)P(k)F  (k)  +  Q(k)  . 

The  Kalman  filter  algorithm  is  implemented  by  recursive¬ 
ly  computing  equations  (2.16),  (2.19),  (2.20),  (2.21),  and 


(2.22). 


Figure  2.3  is  a  diagram  of  the  Kalman  filter.  The 


algorithm  is  initiated  with  the  initial  conditions 


x(0)  =  E[x( 0 ) ] 


and , 


P(0) 


T 

E[ (x(O)-x(O)  ) (x(O)-x(O)  )  ]. 


(2.23a) 


(  2 . 23b) 


In  the  event  that  a  controlling  input  or  a  deterministic 
disturbance  u(k)  is  applied  to  the  system,  the  only  change 
in  the  above  algorithm  is  the  state  model  and  the  a  poster¬ 
iori  estimate.  They  become  [Ref.  5:p.  130] 

x(k)  =  F ( k-1 )x ( k-1 )  +  w(k-l)  +  G ( k- 1 ) u ( k- 1 )  (2.24) 

x(k)  =  F(k-l)x(k-l)+G(k-l)u(k-l)  (2.25) 
+K(k) [z(k)-H(k)F(k-l )x(k-l )  ]  . 

where  G(k)  is  a  (n  x  q)  matrix  and  u(k)  is  a  (q  x  1)  vector 
of  input  signals. 

The  problem  of  estimating  a  desired  signal  based  on 
noisy  data  observations  pertains  to  such  fields  as 
communications,  controls,  and  geophysics.  However,  the 
Kalman  filter  can  also  be  applied  to  these  deconvolution 
problems.  As  an  example,  the  discrete  Kalman  inverse  filter 
will  now  be  applied  to  an  exploration  seismology  problem. 

Typically,  in  the  search  for  underground  oil  and  gas 
deposits  ,  a  vibratory  signal  source  generates  a  pulse  of 
energy  which  is  transmitted  into  the  earth.  In  modern 


seismology,  the  shape  of  this  source  wavelet  can  be 


carefully  controlled;  it  is  chosen  so  that  it  contains  only 
those  frequencies  which  are  transmitted  best  by  the  earth. 
As  the  source  wavelet  propagates  through  the  earth,  it 
encounters  many  different  layers  with  various  acoustic 
impedances.  At  these  layers,  both  partial  reflection  and 
refraction  occur,  creating  numerous  transmission  paths.  The 
received  seismic  signal  at  the  surface  is  composed  of  many 


overlapping  reflected  wavelets.  Therefore,  the  seismic 
trace  can  be  represented  as  the  convolution  of  the  original 
source  wavelet  with  an  impulse  train  representing  the 
various  layers  of  the  earth.  Moreover,  the  seismic  trace  is 
contaminated  by  measurement  noise  and  by  the  phenomena  of 
ghost  reflections  and  reverberations.  [Ref.  6:pp.  14-15] 


The  seismic  trace  can  be  described  mathematically  by 


z(t)  =  s(t,T)*r(t)  +  v(t) 


-r 

j  t. 


s(t,T)r(T)dT  +  v(t)  (2.26) 


where  z(t)  =  measured  seismic  trace 

s(t,T)  =  finite  duration,  time  varying  wavelet 
r(t)  =  reflectivity  function  of  earth’s  structure 
v(t)  =  measurement  noise. 


The  seismologist  must  extract  the  structure  of  the  earth, 
r(t),  by  analyzing  the  noisy  seismic  data  z(t).  This 
process  of  removing  the  wavelet  shape  from  the  trace  and 
leaving  behind  the  impulse  train  representing  the  reflected 
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wavelet’s  strength  and  arrival  time  is  that  of 
deconvolution.  CRUMP  [Ref.  3:pp.  6-7]  shows  that  if  the 
seismic  signal  is  sampled  at  discrete,  uniformly  spaced 
intervals,  and  if  s(t,T)  and  r(t)  are  assumed  to  be  causal, 
then  z(t)  is  represented  by 
J 

z(k)  =  Y  [ s ( k , k-i  +  1 ) r ( k-i  +  1 )  ]  +  v(k)  (2.27) 

i  =  l 

where  the  sample  number  k  =  1,2,3,...,  and 

L  =  length  of  the  wavelet  given  in  number  of  samples 
J  =  k  for  k<L 
=  L  for  k>_L  . 

When  M  traces  of  K  samples  in  length  are  available  for 
processing  then  z(k)  becomes  a  (M  x  1)  vector  where  the  j-th 
component  is  given  by 
J 

z  (k)  =  £  [H  (k))r( k-i+1)]  +  v  (k)  (2.28) 

j  i  =  1  ji  j 

for  j  =  1,2,...,M  and  k  =  1,2,...,K.  This  assumes  that  the 
reflectivity  function  is  the  same  for  each  trace,  while  the 
shape  of  the  exciting  wavelet  may  vary  from  trace  to  trace. 
The  time-varying  wavelet  values  are  contained  in  the  (M  x  L) 
matrix  H(k)  where  the  j-th  row  contains  the  L  samples  of  the 
wavelet  which  generates  the  j-th  trace: 

H  (k)  =  s  (k, k-i  +  1)  .  (2.29) 

j  i  j 

In  vector  form,  the  equations  become 


z  ( k  )  =  H{ k ) x( k )  +  v  (  k  ) 


(2.30) 


where  the  state,  measurement,  and  noise  vectors  are  given  by 

T 

x(k)  =  [ r ( k ), r ( k- 1 ),..., r ( k-L+ 1 ) ]  (2.31a) 

T 

z(k)  =  [z  (k),z  (k),...,z  ( k )  ]  (2.31b) 

1  2  M 

T 

v ( k )  =  [v  ( k ) , v  ( k ),..., v  (k)]  .  (2.31c) 

12  M 

This  is  the  Kalman  filter  measurement  model.  Now  the  state 
model  must  be  determined. 

The  state  model  is  arrived  at  by  assuming  a  general 
relationship  for  the  reflection  coefficients  of  x(k).  The 
assumed  relationship  is  the  autoregressive  equation 
L 

r(k)  =  V  [b  (k-l)r(k-i)]  +  w(k-l)  .  (2.32) 

i  =  1  i 

Comparing  equation  (2.32)  with  the  state  vector  x(k)  yields 
the  state  model 


x(k)  =  F ( k , k- 1 ) x ( k- 1 )  +  w(k-l)  (2.33) 


where  the  (M  x  L)  transition  matrix  and  the  (M  x  1)  system 
noise  vector  are  given  by 


F ( k , k- 1 ) 


b  ( k- 1 )  b  (k- 1 )  ... 

1  b  (k-1 ) 

1  2 

L 

I 

:  o 

i 

i 

(  2 .34a  ) 


p>\ 


K 

f.v 

|.V 


T 

w(k-l)  =  [w(k-l),0,0f...,0]  , 


(  2  .  34b) 


X  is  the  identity  matrix,  and  0  is  the  null  vector. 
Equations  (2.30)  and  (2.32)  provide  the  measurement  and 
system  models  for  implementing  the  deconvolution  via  the 
Kalman  filter.  CRUMP  discusses  methods  by  which  to  obtain 
numerical  values  for  the  reflection  coefficient  vector  b(k) 
and  the  time-varying  wavelet  sample  matrix  H(k)  [Ref.  3 : pp 
8-11].  Once  these  matrices  are  determined,  the  recursive 
Kalman  filter  removes  the  effects  of  s(t,T)  from  r(t)  and 
generates  the  state  estimate  x(k)  which  provides  L  samples 
of  the  desired  reflectivity  function  at  each  time  k. 

It  is  interesting  to  note  that  both  the  Kalman  and 
Wiener  filters  are  minimum  mean-square  error  estimators, 
both  require  the  same  a  priori  knowledge  of  the  process  to 
be  estimated,  and  that  both  yield  identical  estimates. 
However,  the  Kalman  filter  does  have  distinct  advantages 
over  the  Wiener  filter.  First,  due  to  the  matrix  form  of 
the  state  space  equations,  the  Kalman  filter  has 
multichannel  capability  and  is  equivalent  to  a  bank  of 
optimal  estimators.  Moreover,  the  Kalman  filter  is  ideally 
suited  to  computer  implementation  due  to  its  discrete  and 
recursive  characteristics.  [Ref.2:pp.  268-269] 

The  discrete  least-squares  approach  which  foil ows  is  a 
viable  alternative  to  the  Kalman  filter  algorithm, 
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particularly  when  state  space  model  equations  are  not 
available  or  applicable. 

C.  LEAST-SQUARES  INVERSE  FILTER 

When  a  linear  system,  H(z),  is  excited  by  the  an  input 

signal  x(n) ,  the  output  y(n)  is  defined  by  the  convolution 

relationship  y(n)  =  x(n)  *  h(n).  Deconvolution  involves 

finding  the  inverse  filter  G(z)  such  that  H(z)G(z)  =  1. 

(Note  that  the  discrete  time  domain  is  related  to  the  fre- 

jwT 

quency  domain  through  the  equation  z  =  e  ,  where  T  is 
the  discrete  sampling  interval.)  This  condition  transforms 
to  h(n)  *  g(n)  =  d(n)  in  the  discrete  time  domain;  h(n)  and 
g(n)  are  the  impulse  responses  of  the  filters  H(z)  and  G(z), 
respectively,  and  d(n)  is  the  unit  impulse  function.  If 

this  condition  is  met,  then  the  original  input  signal  x(n) 
is  recovered  at  the  output  of  the  inverse  filter. 

The  transfer  function  of  the  causal  system  is  defined  by 
the  infinite  series  obtained  by  taking  the  one-sided  z- 
transform  of  the  system’s  impulse  response: 

co  -i 


H(z)  =  Y(z)/X(z)  =  £  h(i)z 

i  =  0 


(2.35) 


H(z)  can  be  determined  by  inserting  a  known  sequence  x(n) 
into  the  system,  measuring  the  output  sequence  y(n),  and 
then  manipulating  their  z-transforras .  The  inverse  filter 
G(z)  is  then  computed  by  carrying  out  the  polynomial 
division 
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G  (  z  ) 


1/H(z) 


-1  -2 

l/[h(0)+h( 1 )z  +  h ( 2  )  z  +...] 


(2.36) 


-1  -2  -M 

=  g ( 0 )  +  g(l)z  +  g ( 2 ) z  +  ...  +  g ( M ) z  +  ... 

and  truncating  to  M+l  terms  if  necessary.  If  the  exact 

inverse  filter  is  approximated  by  truncating  G(z)  to  order 

M,  then  its  impulse  reponse  is  given  by  the  sequence  g(n), 

n=0 , 1 , 2 , . . . , M .  Furthermore,  if  the  original  system’s 
impulse  response  is  represented  by  h(n),  n=0 , 1 , 2 , . . . N ,  then 

M 

d ( n )  =  g ( n )  *  h ( n )  =  V  g(m)h(n-M)  (2.37) 

m  =  0 

for  0  n  <_  N+M.  This  approximation  to  the  impulse  function 
improves  as  the  order  M  of  the  inverse  filter  is  increased. 

Now  the  stability  of  the  inverse  filter  will  be 
addressed.  If  H(z)  has  all  its  zeros  inside  the  unit  circle 
in  the  complex  z-plane,  it  is  referred  to  as  a  minimum-delay 
polynomial;  the  corresponding  sequence  h(n)  is  called  a 
minimum  phase-lag  sequence.  This  is  a  sufficient  condition 
to  guarantee  that  H(z)  has  a  stable  inverse,  because  the 
zeros  of  H(z)  become  the  poles  of  G(z),  and  G(z)  is  a  stable 
filter  if  all  its  poles  lie  within  the  unit  circle. 

Maximum-  and  mixed-delay  signals  are  obtained  by 

* 

transforming  the  zeros  of  H(z)  from  z  to  1/z  where  the 

i  i 

superscript  represents  the  complex  conjugate  operation. 

A  maximum-delay  sequence  has  all  of  its  zeros  outside  the 
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unit  circle,  while  the  mixed-delay  sequence  has  zeros  inside 
and  outside  the  unit  circle. 

If  H(z)  has  N  zeros,  then  transforming  these  zeros 

N 

results  in  at  most  2  distinct  sequences.  Of  note,  each  of 

these  minimum,  maximum,  and  mixed-delay  signals  have  the 

2  * 

same  magnitude  spectrum:  |H(w)|  =  H(w)H  (w)  [Ref.  l:p.  98]. 

However,  they  do  have  distinct  phase  spectra  [Ref.  4:p. 

175].  The  maximum-delay  polynomial  can  be  written  as 

R  *  *  -1  *  -N 

H  (z)  =  h  (N)  +  h  (N-l)z  +...+  h  <0)z  (2.38) 

R 

The  so  called  reverse  polynomial  H  (z)  is  a  conjugated, 

reflected,  and  shifted  version  of  H(z).  The  corresponding 

R  *  *  * 

maximum  phase-lag  sequence  is  h  =  [h  (N),h  (N-l ),..., h  ( C ) } 
[Ref.  6:p.  72]  While  the  minimum-delay  filter  has  a  causal, 
stable  inverse  consisting  only  of  a  memory  function,  the 
maximum-delay  filter  has  an  inverse  which  consists  only  of  a 
stable,  noncausal ,  anticipation  function.  The  stable 
inverse  of  a  mixed-delay  function  consists  of  both  memory 
and  anticipation  functions.  Filters  with  nonvanishing 
anticipation  components  are  noncausal;  they  cannot  work  in 
real  time  since  the  future  values  of  the  filter  input  are 
not  available  for  processing.  This  problem  can  be 
circumvented  if  the  entire  signal  is  first  recorded  prior  to 
analysis;  then  the  required  future  input  data  is  available. 


[Ref.  4 : p .  87] 


The  energy  distribution  in  minimum, 


mixed,  and  maximum 


delay  signals  will  now  be  examined.  Since  each  of  these 
signals  have  an  identical  magnitude  spectrum,  they  also 
have  the  same  total  energy.  However,  although  the  total 
energy  is  the  same,  the  rate  at  which  the  energy  builds  up 
differs  for  the  various  sequences.  Parseval’s  theorem 
states  that  the  total  energy  in  a  signal  is  given  by 


tr 

I  H  ( w  ) 


2 

dw/ ( 2n ) 


M  2 

£  l  h  (  m  )  | 
m  =  0 


(2.39) 


If 


the  partial 
n 


P(  n) 


energy 

2 

I  h  ( m )  | 


is  defined  as 


(2.40) 


then  it  can  be  shown  that  the  energy  builds  up  quickest  in 
the  minimum-delay  sequence,  and  that  it  builds  up  the 
slowest  in  the  maximum-delay  sequence  [Ref.  6:pp.  75-76].  In 
other  words,  the  minimum-delay  signal  makes  its  impact  as 
soon  as  possible  since  its  energy  is  concentrated  at  the 
front  of  the  sequence.  The  maximum-delay  signal  makes  its 
major  impact  at  a  later  time  since  its  energy  is 
concentrated  at  the  end  of  the  sequence.  The  energy  curves 
associated  with  all  the  possible  mixed-delay  signals  lie 
between  these  two  extremes.  Finally,  it  can  be  shown  that 
the  convolution  of  two  minimum-delay  sequences  with  one 
another  results  in  a  minimum-delay  sequence.  The  convolution 
of  maximum-delay  signals  results  in  a  maximum-delay 


sequence.  Convolution  involving  any  other  combination  of 
sequences  results  in  a  mixed-delay  sequence.  Of  note,  the 
resulting  maximum  and  minimum-delay  sequences  are  the 
reverse  of  each  other.  [Ref.  6:pp.  73-74] 

The  method  described  above  of  finding  an  approximate, 
finite  length  inverse  filter  consisted  of  simply  truncating 
the  exact  inverse  filter  found  by  polynomial  division.  It 
was  seen  that  the  inverse  filter  G(z)  attempted  to  transform 
the  impulse  response  of  H(z)  into  a  unit  impulse  located  at 
the  origin.  This  can  be  thought  of  as  an  attempt  by  G(z)  to 
undo  the  blurring  effect  of  H(z)  (i.e.,  H(z)  "blurs"  the 
impulse  d(n)  into  the  impulse  response  h(n)).  If  the  input 
to  H(z)  is  designated  x{n),  and  if  the  output  of  G(z)  is 
x(n),  then  the  error  of  the  approximated  inverse  filter  is 

e  ( n )  =  x  ( n )  -  x(n)  for  0  <_  n  <  N+M.  (2.41) 

The  error  energy  is  defined  by 
N+M  2 

J  =  £  e  (n)  (2.42) 

n  =  0 

For  the  polynomial  division  /  truncation  method,  J  decreases 
as  the  order  M  of  G(z)  increases  [Ref.  6:p.  136],  Seeking  to 
minimize  the  error  energy  J  with  respect  to  the  inverse 
filter  coefficients  leads  to  another  approach  for  finding  an 


approximate  inverse  filter. 


Filters  which  minimize  the  mean-square  error  J  are 


called  least  error  energy  or  least-squares  inverses.  In 


general,  the  desired  output  sequence  x(n)  of  the  inverse 


filter  G ( z )  could  be  of  any  shape.  G(z)  is  then  called  a 


waveshaping  filter.  When  applied  to  deconvolution,  the 


desired  output  of  the  inverse  filter  is  a  unit  impulse 


d(n-i),  where  i  =  0 , 1 , 2 , . .  .  ,  N+M  defines  the  lag  of  the  digi¬ 


tal  inverse  filter.  In  this  case,  G(z)  is  called  an  i-th 


delay  spiking  filter  since  it  tries  to  condense  the  system 


impulse  response  h(n)  into  a  spike  with  i  delays.  Since  i  = 


0 , 1 , 2 , . . . , N+M ,  there  are  N+M+l  possible  spiking  filters.  As 


will  be  discussed,  there  are  preferred  values  of  the  lag  i 


for  minimizing  J,  depending  on  the  phase  characteristics  of 


h(n).  The  spiking  filter  problem  is  shown  in  Figure  2.4. 


It  is  convenient  to  restate  the  deconvolution  problem  in 


a  matrix  form  based  on  the  Yule-Walker,  or  autocorrelation, 


method  [Ref.  l:pp.  243-245].  This  is  accomplished  by 


defining  the  (M+N+l  x  M+l)  system  output  matrix  Y,  the 


(M+l  x  1)  impulse  response  vector  &,  and  the  (N+M+l  x  1) 


input  vector  x  as  follows: 


< - M  ZEROS 

y  (  0  )  0  0 

y( l)  y(0)  0 

y(2)  y ( 1 )  y(0) 


y(N)  y(N-l)  y ( N-2 ) 
0  y(N)  y(N-l) 


0 

y(N-M) 


(2.43a) 


y(N) 

0 


y  ( N ) 


=  [g(0) ,  g<  1  )  , 


(2.43b) 


=  [x(0),  x(l),  ...  ,x(N+M)] 


( 2 .43c ) 


The  above  notation  is  for  the  general  case  of  the 


waveshaping  filter.  For  the  special  case  of  the  spiking 


filter,  the  y(n)  components  of  the  system  output  matrix  Y 


are  replaced  by  the  impulse  response  h(n).  Also,  the 


desired  signal  x  reduces  to  an  impulse  d(n-i). 


equations  (2.37),  (2.41),  and  (2.42)  can  be  rewritten  as 


A 

x  =  Yg 


(2.44a) 


=  x  -  x  ,  and 


(2.44b) 


( 2 . 44c  ) 


As  discussed  in  the  introduction  to  this  chapter,  the 


1  eas t - squa res  criterion  is  satisfied  by  minimizing  J  with 


respect  to  the  inverse  filter  coefficients  G.  This  results 


in  the  normal  equations 
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T  T 

Y  Yg_  =  Y  x 


(2.45) 


which  can  be  solved  for  g..  By  recognizing  that  Y  Y  is 

equivalent  to  the  (M+l  x  M+l)  sampled  autocorrelation  matrix 
T 

R  =  EtYZ  1  where  the  M+l  length  .data  vector  is  y  =  [  y(0), 

T  T 

y  (  1  y( N ), 0 , 0 ,... 0 ]  ,  and  that  Y  x  is  a  length  (M+l) 

cross-correlation  column  vector  r,  it  can  be  seen  that 


T  -1  T  -1 

g  =  (Y  Y)  Y  x  =  R  r 


(2.46) 


where  R  can  be  evaluated  efficiently  by  Levinson’s 
algorithm.  The  actual  filter  output  is  then 


-1  T 

x  =  Yg  =  (YR  Y  )x  =  Px 


(2.47) 


T  -1  T 

where  P  =  Y(Y  Y)  Y  is  a  square  (N+M+l)  dimensioned  matrix 
called  the  projection  or  performance  matrix.  The  filter's 
performance  improves  as  P  approaches  the  identity  matrix. 
Now,  the  estimation  error  and  cost  function  are  written  as 


=  X  -  X  =  x(I  -  P) 
T  T 

=ee=x(I-P)x 


(2.48) 


(2.49 


Since  in  the  deconvolution  problem  x  is  the  spike  with  i 

delays,  the  inverse  filter  output  x  is  actually  the  i-th 

-1  T 

column  of  the  P  matrix.  Also,  since  £  =  R  Y  ,  the 

coefficients  of  the  i-th  spiking  filter  are  contained  in  the 

-1  T 

i-th  column  of  the  matrix  R  Y  .  Moreover,  the  energy 


1 


error  of  the  i-th  spiking  filter  reduces  to  J  =  1  -  P(i,i). 
Therefore,  in  order  to  realize  the  best  inverse  function 
(i.e.,  minimize  J),  select  the  delay  i  for  which  P(i,i)  is 
largest.  For  a  chosen  lag  i,  the  filter’s  performance  also 
improves  as  the  order  M  of  the  inverse  filter  is  increased. 

Now  let  J ( i )  represent  the  estimation  error  for  the  i-th 
spiking  filter.  The  grand  sum  of  squared  errors  is  defined 
as  V  =  J(0)  +  J(l)  +...+  J(M+N).  It  can  be  shown  that  V  = 
(M+N+l)  -  ( M+ 1 )  =  N,  where  N  is  the  order  of  the  system  H(z) 
[Ref.  4:p.  198].  Therefore,  V  is  independent  of  order  M. 

For  sufficiently  long  spiking  filters,  the  optimal  value 
of  the  delay  i  depends  on  the  phase  characteristics  of  the 
signal  h(n)  and  the  choice  of  the  lag  is  governed  by  the 
following  rules.  If  h(n)  is  a  minimum-delay  input,  the 
spiking  filter  should  have  zero  delay,  i  =  0.  This  says 
that  for  a  signal  with  its  energy  concentrated  towards  the 
front  of  the  sequence,  it  is  easiest  to  condense  it  to  a 
unit  impulse  at  the  origin.  If  the  signal  is  a  maximum- 
delay  input,  the  maximum-delay  spike  i  =  N+M+l  should  be 
selected.  If  h(n)  is  a  mixed-delay  signal,  then  the  i 
corresponding  to  the  largest  P(i,i)  should  be  chosen.  [Ref 
6 . : p .  152] 

Under  certain  conditions,  the  estimation  error  will  go 
to  zero  as  the  length  of  the  inverse  filter  tends  to 
infinity.  First,  as  previously  discussed,  if  h(n)  is  a 


minimum-delay  sequence  then  an  exact  inverse  can  be  found 
through  polynomial  division.  The  zero-delay  spiking  filter 
approaches  this  exact  inverse  filter  G(z)  as  the  number  of 
terms  M+l  increases.  Therefore,  the  estimation  error  goes' 
to  zero  as  M  goes  to  infinity.  The  second  condition  is,  if 
h(n)  is  not  a  minimum-delay  sequence,  J  will  approach  zero 
as  1/M  if  the  lag  i  of  the  spiking  filter  is  chosen  to  be 
sufficiently  large.  [Ref.4:pp.  200-201] 

Up  until  now,  the  effects  of  measurement  noise  and 
imperfect  knowledge  of  the  distorting  function  H(z)  have 
been  neglected.  If  noise  is  introduced,  then  the  output 
y(n)  of  the  system  H(z)  driven  by  input  signal  x(n)  becomes 

y ( n )  =  h(n)*x(n)  +  v(n)  (2.50) 

where  v(n)  is  assumed  to  be  zero-mean  white  noise  with 
variance  Q.  If  this  signal  is  passed  through  the  previously 
determined  inverse  filter  G(z),  the  filter  output  becomes 

x(n)  =  g(n)*y(n)  =  g ( n ) *h ( n  )  *x ( n )  +  g(n)*v(n)  (2.51) 

=  d(n)*x(n)  +  u(n)  =  x(n)  +  u(n) 

where  u(n)  =  g(n)*v(n)  is  the  filtered  noise  signal  and  d(n) 
is  the  unit  impulse  function.  The  variance  of  u(n)  is: 
2  M2  T 

E[u  ( n )  ]  =  Q  £  g  (n)  =  Q&  g_  (2.52) 

n  =  0 

[Ref.  l:p.  248].  This  variance  may  be  larger  than  the 


original  noise  variance  Q.  For  example,  if  h(n)  is  a  low 
frequencey  signal,  then  g(n)  must  be  very  spiky  in  order  to 
compress  h(n)  into  a  high  frequency  content  spike.  As  a 
result,  g<n)  will  likely  have  values  greater  than  one. 
Therefore,  the  variance  of  u(n)  will  be  larger  than  Q.  This 
further  degrades  the  estimate  x(n). 

To  compensate  for  this  filtered  noise,  the  minimization 
criterion  is  modified  so  that 

N+M  2  T 

J  =  Y  (d(n)  -  g(n)*h(n))  +  XQg_  g_  (2.53) 

n  =  0 

where  A  is  a  positive  parameter.  The  first  term  of  the  cost 
function  tries  to  produce  a  good  inverse  filter  whereas  the 
second  terra  tries  to  reduce  the  output  noise.  If  A  is  large, 
noise  reduction  is  emphasized  at  the  expense  of  obtaining  a 
good  inverse  function.  A  small  A  emphasizes  finding  a  good 
inverse  function  g(n),  and  there  is  little  output  noise 
reduction . 

Using  this  new  cost  function,  the  resulting  normal 
equations  are  given  by 
T  T 

(Y  Y  +  XQI )a  =  Y  x  .  (2.54) 

Comparing  this  with  equation  (2.45)  reveals  that  the  main 
diagonal  elements  of  the  sampled  autocorrelation  matrix  R 
have  been  modified  by  an  additive  term:  R(0)  becomes  R(0)  + 
AQ .  If  the  Backus-Gilbert  " prewhi ten ing"  parameter  epsilon 


is  introduced,  where 


e  =  AQ/R(0)  ,  (2.55) 

then  the  diagonal  elements  are  written  as  {1  +£  )R(0).  Even 

very  small  values  for  the  BG  parameter  have  a  beneficial 

T 

effect  of  stabilizing  the  inverse  of  the  matrix  (Y  Y  +  AQD* 
[Ref.  1 : p .  246] 

D.  LATTICE  FILTER 

1 •  Introduct ion 

The  lattice  filter  is  another  optimal  least-squares 
predictor  which  can  be  applied  to  the  linear  deconvolution 
problem.  The  lattice,  or  ladder  filter  derives  its  name 
from  the  cascade  form  of  its  signal  flow  graph.  Basically, 
the  lattice  filter  provides  an  alternative  to  the 
transversal  filter  for  modeling  a  signal.  It  can  be  viewed 
as  a  Gram-Schmidt  orthogonalization  of  the  incoming  data. 
This  will  be  elaborated  upon  in  subsequent  paragraphs.  To 
date,  much  of  the  work  done  with  lattice  filters  has  been  in 
the  areas  of  speech  analysis/synthesis,  seismology,  and  in 
high-resolution  spectral  estimation  [Ref.  7:p.  841].  Prior 
to  developing  the  lattice  filter  equations,  key  mathematical 
concepts  which  apply  to  this  discussion  will  be  reviewed. 
2 .  Mathematical  Background 

A  central  concept  behind  the  development  of  the 
lattice  filter  is  that  of  orthogonality.  Two  vectors  are 


orthogonal  if  their  inner  product  is  zero.  The  geometric 
interpretation  of  this  is  that  the  vectors  are  at  right 
angles  to  one  another.  As  an  example,  if  the  random  varia¬ 
bles  u  and  v  are  the  basis  vectors  of  a  two-dimensional 

space,  they  are  orthogonal  if  and  only  if  their  inner 

product  <u,v>  =  E[uv]  =  0.  In  the  case  where  the  linear 
space  is  spanned  by  the  random  variables  of  a  data  vector, 
two  vectors  are  orthogonal  if  the  corresponding  random 

variables  are  uncorrelated  and  if  one  or  both  have  zero 
mean  [Ref.  8:p.  92]. 

A  related  theorem  is  the  orthogonal  decomposition 
theorem,  which  states  that  any  random  variable  may  be 
decomposed  uniquely  with  respect  to  a  subspace  S  into  two 
mutually  orthogonal  parts,  one  part  which  is  parallel  to  S 
(i.e.,  lies  in  S)  and  the  other  part  orthogonal  to  S.  That 

A  A 

is,  £  can  be  written  y_  =  Z  +  §.  where  e  is  orthogonal  to  y_ 

and  to  the  basis  vectors  which  span  the  subspace  S,  and 

where  £  is  defined  by  a  linear  combination  of  the  basis 

vectors.  If  S  is  spanned  by  the  basis  vectors 

M 

A  r— » 

{ u ( 1 ), u ( 2 ),..., u ( M )} ,  then  y  =  )  a(i)u(i)  where  it  can  be 

i  =  1 

shown  that  the  coefficients  are  given  by  the  equation 

2  -1 

a ( i )  =  E[yu(i)]E[u  (i)]  .  [Ref.  l:p.  13] 

The  orthogonal  projection  theorem  adds  to  this  by 
stating:  the  orthogonal  projection  £  of  a  vector  onto  a 

linear  subspace  S  is  that  vector  in  S  which  lies  closest  to 
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Y_  with  respect  to  the  distance  induced  by  the  inner  product 
of  the  vector  space  [Ref.  l:p.  15].  Restated,  this  simply 
means  that  £  is  the  best  estimate  of  y;  that  can  be  made  by  a 
linear  combination  of  the  basis  vectors  of  S  in  a  minimum 
mean-squared  error  sense.  Figure  2.5  shows  the  projection 
of  v  into  the  subspace  S. 

Another  important  concept  in  developing  the  lattice 
filter  is  that  of  Gram-Schmidt  orthogonalization .  The  Gram- 
Schmidt  procedure  is  a  recursive  orthogonalization  process 
which  generates  a  set  of  mutually  orthogonal  basis  vectors 
{u ( 1 ), u ( 2 ),..., u ( M ) }  from  a  given  set  of  basis  vectors 
{ y ( 1 )i y ( 2 ),..., y ( M )} .  The  procedure  is  initialized  by 
letting  u( 1 )  =  y(l).  Next,  y(2)  is  decomposed  with  respect 
to  u(l).  That  part  of  y(2)  which  is  orthogonal  to  u(l) 
becomes  u(2).  Next,  y(3)  is  decomposed  with  respect  to  the 
subspace  spanned  by  {u(l),u(2)}.  That  part  of  y(3)  which  is 
orthogonal  to  this  subspace  becomes  u(3).  In  general,  the 
new  set  of  orthogonal  basis  vectors  is  defined  by 

n- 1 

u(n)  =  y(n)  -  £  b(n,i)u(i)  (2.56) 

i  =  1 

2  -1 

for  2  n  <.  M  and  where  b(n,i)  =  E  [  y  (  n  )  u  (  i  )  ]  E  [  u  (i)] 

Using  the  orthogonal  decomposition  theorem,  this  is 
equivalent  to  y(n)  =  y(n)  +  j(n)  where  y(n)  is  the  best 
estimate  of  the  n-th  component  of  the  data  vector  based  on 
the  previous  n-1  components  and  u(n)  represents 


the 


estimation  error.  If  b(n,n)  =  1  is  introduced,  then 


This 


y  (n) 

can  be 


n 

b(n, i )u( i ) 

1 

written  in  matrix 


for  1  <_  n  <_  M. 

form  as 


(2.57a) 


=  Bu  where  (2.57b) 

T 

Z  =  f y( 1 ) »  y( 2 ) , . . . ,  y ( M) ]  , 

T 

u  =  [u( 1 ) ,  u ( 2  )  ,  . . . ,  u ( M ) ]  , 


B 


1  0  0 

b<  2 ,  1  )  1  0 

b( 3 ,  1  )  b ( 3 , 2  )  1 


0 

0 

0 


b  ( M ,  1  )  b  (  M  ,  2  ) 


This  is  a  convenient  notation  for  representing  the 

transformation  from  a  set  of  correlated  basis  vectors  z  to 

an  uncorrelated  set  of  basis  vectors  u.  The  bases  z  and  y 

span  the  same  M-diraens ional  subspace  S,  but  there  are  no 

redundant  correlations  between  the  basis  vectors  of  set  u. 

Since  the  basis  of  u  is  uncorrelated,  its  components  u(i) 

(i  =  1,2,...,M)  are  referred  to  as  innovations  because  each 

additional  component  contributes  completely  new  information 

to  the  estimate  of  z*  [Ref.  l:pp.  16-18] 

Finally,  the  transformation  z  =  By.  corresponds  to  a 

LU  (lower  upper ) -Cholesky  factorization  of  the  correlation 

T 

matrix  of  z*  By  definition  R  =  E[zz  )•  Substituting  for 

yy 
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Y_  results  in 


K 


T  T  T 

R  =  BE[uu  ]B  =  BR  B 
yy  uu 


(2.58) 


where  B  is  lower  triangular,  B  upper  triangular,  and  R  is 

uu 

a  diagonal  matrix  since  the  basis  vectors  u(i)  are 

uncorrelated  with  one  another. 

3 .  Derivation  of  Lattice  Filter  Equations 

Now  the  lattice  filter  order  update  equations  will 

be  developed.  A  random  signal  y(n)  can  be  modeled  as  the 

output  of  a  causal,  stable,  linear  filter  H(z)  which  is 

driven  by  a  stationary,  uncorrelated,  white  noise  sequence 

{ e ( 1 ) , e ( 2 ) . e ( P ) }  [Ref.  l:p.  30].  A  P-th  order 

autoregressive  model  is  defined  by  H(z)  =  1/A  (z)  where 

P 

-1  -2  -P 

A  (z)  =  l+a(l)z  +a(2)z  +...+a(P)z  .  (2.59) 

P 

The  filter  A  (z)  has  several  names:  prediction  error  filter, 
P 

inverse  filter,  or  analysis  filter.  The  signal  y(n)  can  be 
written  as 


y ( n )  =  e ( n )  -  Y  a(i)y(n-i) 
i  =  1 


(2.60) 


^  i  ^  m  ■  M'P  1 1  y  1 7 w'l  i 


v( n-P+l ) , . . . ,  y(n-l)}  and  e(n)  is  equivalent  to  the  estima¬ 

tion  error.  The  estimate  y(n)  can  be  thought  of  as  the 
projection  of  the  ( P+ 1 ) -dimensional  y(n)  vector  onto  the  P- 

dimensional  subspace  spanned  by  the  components  of  the  data 

T 

vector  Y  =  [ y ( n- 1 ), y ( n- 2 y ( n-P ) ]  .  In  general,  the 
past  values  of  y(n)  are  correlated  with  one  another.  There¬ 
fore,  the  random  variable  components  of  Y  do  not  generally 
form  a  set  of  mutually  orthogonal  basis  vectors.  The  Gram- 
Schmidt  or thogonal i zation  procedure  removes  these 
correlations  and  transforms  Y  into  a  set  of  mutually 
orthogonal  basis  vectors  which  span  the  same  subspace. 
Additionally,  the  predictor  coefficients  a(i)  are  replaced 

by  the  reflection  coefficients  K  .  This  results  in  the 

i 


lattice  filter. 

In  order  to  obtain  the  optimal  estimator,  the 

predictor  coefficients  which  minimize  the  mean-squared 

2 

prediction  error  E(e  ( n  > ]  must  be  determined.  This  problem 
was  discussed  in  the  introduction  to  this  chapter.  The 
resulting  matrix  form  of  the  normal  equations  is 


R  (  0  )  R  (  1  )  R  {  2  )  ...  R  (  P  ) 

1 

Ep 

R  (  1  )  R  (  0  )  R  (  1  )  ... 

a  (  1  ) 

0 

R  (  2  )  R  (  1  )  R  (  0  )  ... 

•  •  •  • 

a  (  2  ) 

0 

• 

•  •  •  • 

•  •  •  • 

R  (  P  )  R  (  0  ) 

• 

• 

a(  P) 

» 

0 

where  R(k)  =  E [ y ( n+k ) y ( n ) ]  and  E  is  the  minimized  mean- 

P 

squared  prediction  error  for  the  P-th  order  filter.  The 
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sample  autocorrelation  matrix  components  are  calculated  from 
the  length  N  data  sequence  by 


N-l-k 

R ( k )  =  (1/N)  V  y ( n+k ) y ( n )  (2.63) 

n  =  0 


for  0  k  P  and  where  P  ±  N-l.  [Ref.  l:p.  150] 

Now  there  are  P+1  equations  and  P+1  unknown  model 

parameters  { a( 1 ) , a( 2 ) , . . . , a ( P ) ; E  }.  The  P+1  equations  can  be 

P 

solved  by  inverting  the  autocorrelation  matrix.  This 
3  2 

requires  0(P  )  operations  and  0(P  )  storage  locations.  The 

P+1  equations  can  be  solved  more  efficiently  by  taking 

advantage  of  the  matrix’s  Toeplitz  structure  by  using 

Levinson’s  algorithm.  Levinson’s  algorithm  reduces  the 

2 

required  number  of  operations  and  storage  locations  to  0(P  ) 
and  O(P),  respectively  [Ref.  l:pp.  150-151],  Being  a 
recursive  procedure,  Levinson’s  algorithm  permits  the 
calculation  of  the  (P+l)-st  order  model  parameters  by  using 
the  previously  determined  P-th  order  model  parameters.  The 
matrix  form  of  the  algorithm  is  given  by 


*  v  ■>  V* 


'S* 


$ 

IV 


a  (  1  ) 
P+1 

a  (  2  ) 
P+1 


a  (P) 
P+1 

a  (P+1) 
P+1 


a  (1) 
P 

a  (2) 
P 


a  (P) 
P 


a  (P) 

P 

a  (P-1) 
P 


a  (  1  ) 
P 


(2.64) 


where 


(2.65) 


Y  a  ( i ) R ( P+ 1 - i ) 
i  =  0  P 

i  =  - 

P+1  P 

Y  a  ( i )R( i ) 

i  =  0  P 


The  P  subscript  on  the  "a"  parameters  specifies  the  P-th 

order  prediction  error  filter  A  (z)  whereas  the  index 

P 

identifies  the  appropriate  term  in  the  A  (z)  polynomial. 

P 

K  is  called  the  (P+l)-st  order  reflection  or  PARCOR 

P+1 

(partial  correlation)  coefficient.  The  PARCOR  coefficient 

K  represents  the  true,  or  direct,  correlation  between 

P+1 

y(n-P-l)  and  y(n)  with  the  effects  of  the  intermediate 
variables  (i.e.,  y ( n-P ), y ( n-P+ 1 ),..., y ( n- 1 ) )  removed.  The 
recursive  form  of  Levinson’s  algorithm  is  written  as 


a  (m)=  a  (m)  -  K  a  (P+l-m)  for  l<m<P,  (2.66a) 
P+1  P  P+1  P 


one-to-one  correspondence  between  the  PARCOR  coefficients 


and  the  coefficients  of  the  transfer  function  A  (z) .  The 

P 

transfer  function  or,  equivalently,  the  autocorrelation 

T 

matrix  R  =  E[y(n)y  (n)]  uniquely  determine  the  reflection 
coefficients  [Ref.  7:p.  829].  Taking  the  z-transform  of 
equation  (2.66)  results  in 


where 


A  (z)  =  A  (z)  -  K  B  (z) 
P+1  P  P+1  P 


-1  -2  -P 

A  (z)  =l+a  (l)z  +a  (2)z  + . . . +a  (P)z 

P  P  P  P 


( 2.67a) 


(2.67b) 


-P  -1 

B(z)=zA(z),  (2.67c) 

P  P 

-1  -2  -(P-1)  -P 

=  a  (P)+a  (P-1 ) z  +a  (P-2)z  +...+a  (l)z  +z 

P  P  P  P 

Due  to  the  effective  folding  about  the  axis  and  the  shifting 

to  the  right  of  the  sequence  A  (z),  the  polynomial  B  (z)  is 

P  P 

called  the  reverse  of  A  (z).  Taking  the  reverse  of  both 

P 

sides  of  equation  (2.67)  yields 


B  (z)  =  z  B  (z)  -  K  A  (z) 
P+1  P  P+1  P 


(2.68) 
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Combining  equations  (2.67)  and  (2.68)  into  matrix  form  gives 


-i 

A  (z) 

1  -K  z 

A  (z) 

P+1 

P+1 

P 

-1 

B  (z) 

-K  z 

B  (z) 

P+1 

P+1 

P 

The  forward  prediction  error  associated  with 

predicting  y(n)  from  the  previous  P  samples  { y ( n-P ) , . . . y ( n- 

1))  is  written  as  e  (n)  =  y(n)  -  y  (n)  where  the  subscript  P 

P  P 

represents  the  filter  order.  In  z-domain  notation,  this 
becomes 


E  (z)  =  A  (z)Y(z)  (2.70) 

P  P 

Now,  the  backward  prediction  error  r  (n)  is  introduced.  It 

P 

is  defined  as  the  error  in  predicting  (  or  actually 
smoothing)  y(n-P)  from  the  future  P  samples  { y ( n-P+ 1 , . . . , 
y(n)}.  It  is  written  as 


r  (n)  =  y(n-P)+a  ( 1 ) y ( n-P+ 1 ) + . . . +a  (P)y(n)  (2.71a) 

P  P  P 

or  in  z-domain  notation  as 


E  (z)  =  B  (z)Y(z)  .  (2.71b) 

P  P 

The  optimal  backward  predictor  coefficients  minimize  the 

2 

mean-squared  smoothing  error  E[r  (n)].  Since  the  optimal 

P 

backward  and  forward  predictor  coefficients  are  mirror 


images  of  each  other,  then  E[r  (n)]  =  E[e  (n)].  That  is, 
the  backward  and  forward  prediction  error  vectors  have  the 
same  norms  [Ref.  8:p.  101]. 

As  will  be  shown,  at  each  instant  n,  the  backward 
prediction  errors  are  mutually  orthogonal  (i.e.,  uncorre¬ 
lated  if  assumed  to  be  zero  mean)  [Ref.  l:p.  170]. 
Therefore,  it  is  the  backward  prediction  errors, 

r  ( n ) , p=0 , 1 , . . . , M ,  where  M  is  the  filter  order,  which  form 
P 

the  new  set  of  basis  vectors.  To  demonstrate  that  the 
backward  prediction  errors  are  mutually  orthogonal,  let  M=3 
and  write  the  corresponding  equations  for  P  =  0, 1,2,3  in 
matrix  form: 


r  (n ) 

0 

1  0 

0 

0 

y(n) 

r  (  n ) 

1 

a  (  1  )  1 

1 

0 

0 

y ( n-1 ) 

r  ( n ) 

2 

a  (2)  a  (1) 

2  2 

1 

0 

y ( n-2 ) 

r  (  n ) 

3 

a  (3)  a  (2) 

3  3 

a  (  1  ) 
3 

0 

y ( n-3 ) 

r 

r(n)  =  Ly(n)  .  (2.72b) 

Note  that  the  first  column  of  L  contains  the  negatives  of 
all  the  reflection  coefficients.  Now  examine  the  covariance 


R 


=  E  [  r  (  n )  r  (n)]  =  E  [  Ly.(  n  )  (  Ly(  n  )  )  ] 


(2.73) 


rr 


T  T  T 

LE[y;(n)y;  (n)]L  =  LR  L 

yy 


Since  the  normal  equations  are  satisfied  within  the  matrix 

products  above,  R  reduces  to  a  diagonal  matrix  which 

rr 

verifies  that  the  components  of  r(n)  are  uncorrelated.  That 
is 


T 

R  =  LR  L  =  D  =  diag { E  ,E  ,E  ,E  }  (2.74) 

rr  yy  0  12  3 

where  E  is  the  minimum  value  of  the  mean-squared  prediction 
P  2 

error  which  is  given  by  E[e  (n)  •]  .  It  can  be  shown  that 
2  P+1 

E  =  ( 1-K  )E  where  the  recursion  is  initialized  with 
P+1  P+1  P  2 

the  value  given  by  E  =  R  (0)  =  E[y  (n)]  [Ref.  8:p.  105]. 

0  yy 

Therefore,  the  prediction  error  decreases  by  a  factor  of 
2 

(1-K  )  from  one  lattice  stage  to  the  next.  Now ,  since  the 

P+1 

elements  of  r(n)  are  uncorrelated,  equation  (2.72)  is 

equivalent  to  the  Gram-Schmidt  orthogonal i zation  of  the  data 

vector  ;y(n).  Also  note  that  rewriting  equation  (2.74)  as 
-1  -T 

R  =  L  DL  corresponds  to  a  LU-Cholesky  factorization  of 

yy 

the  covariance  matrix  of  y(n). 

To  complete  the  derivation  of  the  lattice  recursion 
equations,  multiply  both  sides  of  equation  (2.69)  by  Y(z)  to 
get 


T 


-1 

■ 

E  (z) 

1  -K  z 

E  (z) 

P+1 

P+1 

P 

=  -1 

E  (z) 

-K  z 

E  (z) 

P+1 

P+1 

P 

These  equations  are  transformed  into  the  time  domain  to 
obtain  the  final  result: 


e  (n)  =  e  (n)  -  K  r  (n-1) 
P+1  P  P+1  P 

r  (n)  =  r  (n-1)  -  K  e  (n) 
P+1  P  P+1  P 


These  equations,  which  are  recursive  in  both  time  and  order, 
define  the  signal  flow  graph  of  the  analysis  lattice  filter, 
shown  in  Figure  2.6.  For  a  given  time  instant  n,  the  equa¬ 
tions  are  evaluated  recursively  in  order.  The  inputs  into 

the  first  stage  of  the  lattice  filter  are  e  (n)  =  r  (n)  = 

0  0 

y(n).  The  successive  stages  of  the  filter  develc >  the 

successively  higher  order  forward  and  backward  prediction 

errors.  The  output  from  the  final  stage  yields  the  desired 

M-th  order  forward  prediction  error,  while  all  lower  order 

prediction  errors  are  available  at  intermediate  stages. 

Since  the  backward  errors  are  generated  from  the  y  data 

vector,  the  analysis  lattice  filter  actually  implements 

the  Gram-Schmidt  orthogonal i zat i on ;  all  that  is  required  to 

implement  the  lattice  are  the  reflection  factors 

(K  ,K  . K  }  . 

1  2  M 


( 2 . 76a ) 

(2.76b) 


•'  .  V  "  •  V  “j*  V  V.'-'  V.V  h-  w. 


Equation  (2.76)  can  be  manipulated  to  obtain  the 

lattice  form  of  the  synthesis  filter  H(z)  =  1/A  (z): 

M 

e  (n)  =  e  (n)  +  K  r  (n-1)  (2.77a) 

P  P+1  P+1  P 

r  (n)  =  r  (n-1)  -  K  e  (n)  .  (2.77b) 

P+1  P  P+1  P 

Figure  2.7  shows  the  corresponding  signal  flow  graph.  When 

the  input  to  the  synthesis  filter  is  the  forward  prediction 

error  sequence  e  (n),  the  output  is  the  original  sequence 

P 

y(n).  In  order  for  the  synthesis  filter  to  be  stable  and 

causal,  all  M  zeros  of  the  prediction-error  filter  A  (z) 

M 

must  lie  within  the  unit  circle  in  the  complex  z-plane.  A 

necessary  and  sufficient  condition  for  all  the  zeros  to  be 

inside  the  unit  circle  is  that  the  magnitude  of  each  of  the 

reflection  coefficients  {K  ,K  ,...,K  }  be  less  than  one. 

1  2  M 

[Ref.  1  :  pp  .  168-169  ]__ 

The  reflection  factors  can  be  evaluated  by  several 

methods.  The  various  methods  arise  due  to  different 
definitions  of  the  optimality  criterion.  The  criterion 
used  here  of  minimizing  the  mean-squared  prediction  errors 
leads  to  what  MAKHOUL  calls  the  forward  and  backward  methods 
[Ref.  9].  They  are,  respectively: 

2 


K 

=  E[e  (n)r  (n-1)] 

/ 

E  [  r 

(  n-1  )  ] 

( 2 . 78a ) 

P+1  ,e 

P  P 

P 

2 

K 

=  E[e  (n)r  (n-1)] 

/ 

E  [  e 

<  n  )  ] 

( 2 .78b) 

P+1  ,  r 

P  P 

P 
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Assuming  stationarity ,  and  since  E[r  (n)  ]  =  E(e  (n)  ]  as 

P  P 

previously  discussed,  then  the  forward  and  backward 

reflection  coefficients  are  equal.  That  is  K  =  K  = 

P+1  P+ 1 , e 

K  .  The  Schwarz  inequality  implies  that  K  !  <_  1  for 

P+ 1 , r  P 

each  P=1,2,...,M  [Ref.  8:p.  104].  An  alternate  technique  for 

calculating  the  reflection  coefficients  is  the  geometric 

mean  method.  It  was  introduced  by  ITAKURA  and  SAITO  in 

their  work  of  developing  a  digital  filter  structure  for 

time-domain  speech  analysis  [Ref.  10].  The  corresponding 

PARCOR  coefficient  is  given  by 


(2.79) 


E[e  (n)r  ( n— 1 )  ] 

P  P 

2  2  1/2 
(E[e  (n)]E[r  ( n— 1 ) ] } 

P  P 


These  PARCOR  coefficients  are  guarenteed  to  have  magnitudes 

less  than  one  [Ref.  7:p.  840]. 

A  third  method  was  used  by  BURG  in  the  maximum- 

entropy  method  of  spectral  estimation  [Ref.  11].  This  is 

the  harmonic-mean  method.  The  harmonic-mean  method  seeks  to 

minimize  the  sum  of  the  forward  and  backward  prediction 

2  2 

error  variances,  E[e  (n)]  +  E[r  ( n— 1  )  ] .  Minimizing  this 

P  P 

sum  with  respect  to  the  reflection  coefficients  results  in 


(2.80) 


2E[e  (n)r  (n-1 ) ] 

P  P 

K  =  - 

P+1  2  2 

E[e  (n)  +  E[r  (n-1)] 
P  p 


Again,  the  Schwarz  inequality  verifies  that  K  has  magnitude 

P 

less  than  one,  guaranteeing  that  the  synthesis  fil+^r  is 
causal  and  stable  [Ref.  l:p.  189]. 

4 .  The  Generalized  (Analysis)  Lattice  Filter 

The  preceding  development  of  the  analysis  lattice 
filter  equations  assumed  that  the  data  sequence  { y ( n ) } 
represented  a  time  sequence  with  stationary  statistics.  In 
this  section,  a  more  general  linear  prediction  problem 
will  be  considered.  No  special  assumptions  are  made 
concerning  the  data.  The  data  values  need  not  be  delayed 
versions  of  each  other;  they  need  not  even  represent  a  time 
sequence.  This  is  the  approach  taken  by  LENK  in  developing 
the  generalized  order  update  equations  [Ref.  12:p.  85].  The 
resulting  generalized  form  of  the  lattice  filter  makes  it 
suitable  for  multidimensional  and  nonlinear  signal 
processing  applications. 

Definitions  associated  with  the  normalized  form  of 
the  generalized  lattice  filter  will  now  be  introduced. 
First,  the  components  of  the  length  (M+l)  column  vector 
[ y  A  ] ,  representing  a  single  realization  of  the  random 
process  Y,  are  designated  by  y*  ,  A=  0,1,..., M.  The  forward 


prediction  error  in  estimating  the  (n+l)-st  element  from 


the  previous  N  elements  of  the  data  vector  is  written  as 


where  the  normalized  prediction  error  coefficients  are 
defined  as 


N  N  N 

a  ( n+1 )  =  h  ( n+ 1 )  /  | | e  II 
x  A  n+1 


(2.87a) 


N  N 

b  (n-N)  =  h  (n-N)  /  ||r  | | 

A  A  n-N 


(2.87b) 


A  A  n-N 
[Ref.  1 2 : pp .  85-87] 

Using  the  normalized  form  of  the  generalized 
Levinson  algorithm,  LENK  demonstrated  that  equations  (2.87) 
could  be  updated  recursively  through  the  relation 


(2.88) 


N+1 

N 

a  ,  ( n+1 ) 

n 

a  ( n+1 ) 

A 

=  Q(  K  ) 

A 

N+1 

N  +  1 

N 

b  (n-N) 

b.(n-N) 

where  the  partial  correlation  (PARCOR)  coefficient  is 
n  N  N 


n  N  _N 

K  =  E  { e  r  }  , 

N+1  n+1  n-N 


(2.89) 


and  where 


n  1  11 

$(K  )  =  -  N+1  (2.90) 

N+1  n  2  n 

1  -  (K  )  -  K  1 

N+1  N+1 


The  recursion  is  started  for  order  N=0  with  the  initial  pre¬ 
diction  error  values  for  A=0,1,...,M  given  by  the  vectors 


(2.91a) 


0  n+l 

[ ( n+ 1 ) ]  =  [0, . . . ,0,  1 / I  I y  II  ,0 . 0] 

0  n 

[bx(n)j  =  l/||y  II,  0 . 0]  .  (2.91b) 

Multiplying  both  sides  of  equation  (2.90)  by  the  data  vector 
[y  *  ]  yields  the  desired  error  order  update  equations: 


N+l 

_N 

e 

n+  1 

=  0<K  "  ) 

n+  1 

N+l 

N+l 

N 

r 

n-N 

n-N 

(2.92) 


The  corresponding  signal  flow  graph  for  a  single  lattice 
filter  section  is  shown  in  Figure  2.8.  Figure  2.9  depicts  a 
third  order  generalized  lattice  filter.  [Ref.  1 2 : pp .  94-99] 
LENK  then  proved  that  the  reflection  coefficients 
(i.e.,  PARCOR  coefficients)  could  be  calculated  directly 
from  the  data  vector’s  correlation  matrix  by  utilizing  the 
generalized  Schur  algorithm.  In  LENK’s  notation,  the 
components  of  the  correlation  matrix  are  written  as  R^ 
E{yf*  y^  }•  Then  the  reflection  coefficients  are  given  by 

n-N 

a  (n+l) 
n  N 


,rA. 


K 


(2.93) 


N+l 
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n-N 


(  n-N ) 


N 


where 
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(  2 .94a  ) 


<y  (n+l  )  =  E{e  y  }  =  £  a  <n+l)R 

N+l  n+l  p=0 


and 


_  A.  _N  X  M  N  r 

/?  (n-N)  =  E{r  y  }  =  V  b  (n-N)R 

'N+l  n-N  p  =0 


A 


(2.94b) 


Equations  (2.94)  can  be  updated  through  the  recursion 
relation 


A 

A 

Oi  (n+l) 

Oi  (n+l) 

N+l 

n 

N 

=  (K  ) 

x 

N+l 

_A 

fi  (n-N) 

f~N+l 

fi  (n-N) 

^  N 

(2.95) 


To  start  the  recursion  at  order  N=0,  the  values  of  the 

parameters  for  ,\=0,1,...,M  are  given  by  the  vectors 

X  n+l  -1  ( n+ 1 ) 0  ( n+ 1 ) 1  (n+l)M 

[(X  (n+l)]  =  |  |y  II  [R  ,  R  ,  .  .  .  ,R  ]  (2.96a) 

0 

X  n-ln0nl  nM 

[fi  (n)]  =  I  I  y  II  [R  ,R  ,  .  .  .  ,R  ]  .  (2.96b) 

'o 

(Ref.  1 2 : pp .  100-101] 

5 .  Lattice  Filter  Advantages 

The  lattice  filter  has  several  advantages  compared 

to  the  direct,  or  transversal,  form  of  the  prediction  error 

filter  A  (z).  First  of  all,  due  to  the  built-in 
M 

orthogonal ization  incorporated  into  the  lattice,  successive 


lattice  stages  are  decoupled. 


Therefore,  the  reflection 


coefficients  K  ,  P=1,2,...,M  are  independent  of  the  filter’s 
P 

final  order.  The  M-th  order  least  squares  prediction  filter 

contains  all  the  prediction  error  filters  of  lower  order. 

Restated,  the  first  P  sections  of  the  M-th  order  filter  form 

the  P-th  order  prediction  filter.  Therefore,  lattice  stages 

may  be  added  or  subtracted  from  the  existing  lattice  filter 

without  having  to  recalculate  the  already  determined 

reflection  coefficients.  In  contrast,  when  the  order  of  the 

transversal  filter  is  changed,  all  the  A  (z)  prediction 

M 

error  filter  coefficients  must  be  recalculated.  As  a 
consequence,  to  specify  all  filter  orders  up  to  M  requires 
storage  of  M(M+l)/2  predictor  coefficients,  while  only  M 
reflection  coefficients  need  to  be  stored.  [Ref.  7:p.  830] 

Another  advantage  of  the  lattice  over  the 
transversal  filter  is  that  the  output  of  each  lattice  stage 
is  a  least-squares  prediction  error.  Therefore,  if  the 
desired  order  of  the  predictor  is  not  known  in  advance,  the 
output  of  the  various  stages  can  be  monitored  to  determine 
what  filter  order  is  adequate.  [Ref.  8:p.  106] 

Due  to  the  quantization  inherent  in  implementing 
digital  filters,  round  off  noise  is  introduced.  Studies 
have  shown  that  the  performance  of  the  lattice  filter  is 
far  superior  to  that  of  transversl  filters  for  finite  word 
length  computations.  Furthermore,  the  filter  zeros  are  less 
sensitive  to  quantization  errors  in  the  reflection 


coefficients  than  for  quantization  errors  in  the  transversal 
filter  coefficients.  ’  Finally,  since  the  magnitude  of  the 
reflection  coefficients  is  less  than  one,  it  simplifies  the 
task  of  establishing  the  quantizer  overload  point.  [Ref. 
8 : p .  106] 

Another  advantage  the  lattice  filter  holds  over  the 
transveral  filter  is  that  the  lattice  filter  is  minimum 
phase  (i.e.,  all  its  zeros  are  inside  the  unit  circle  so 
that  the  synthesis  model  is  stable)  if  and  only  if  the 
magnitude  of  the  reflection  coefficients  is  less  than  one. 
There  is  no  comparable  test  for  the  transversal  filter 
coefficients  to  determine  whether  the  synthesis  model  is 
stable.  [Ref.  8:p.  106] 

6 .  Linear  Lattice  Filter  Applied  to  Deconvolution 

The  transfer  function  of  the  lattice  filter  is 

determined  by  the  reflection  coefficients.  In  turn,  the 

values  of  these  reflection  coefficients  are  uniquely 

determined  by  the  prediction  error  filter  transfer  function 

A  (z),  or  equivalently  by  the  autocorrelation  sequence  R(k), 
M 

k  =  0,1,.  ...M  [Ref.  7:p.  829].  Similarly,  equations  (2.70) 

and  (2.71b)  indicate  that  the  transfer  functions  from  the 

input  y(n)  to  the  outputs  e  (n)  and  r  (n)  are  A  (z)  and 

M  M  M 

B  (z),  respectively.  Therefore,  the  analysis  lattice  filter 
M 

is  equivalent  to  the  whitening  filter  A  ( z ) ;  that  is,  it  is 

M 

=  1/A  (z).  This  is  the 
M 


the  inverse  of  the  system  model  H(z) 


a 


VM 


a 


basis  for  using  the  lattice  filter  in  a  deconvolution 
application . 

In  order  to  recover  the  input  signal  x(n)  from  the 

system’s  output  signal  y(n),  it  is  necessary  to  determine 

the  inverse  of  the  system  transfer  function  H(z).  The 

initial  step  is  to  conduct  an  "identification  experiment"  to 

estimate  the  system’s  parameters  (either  the  autoregressive 

filter  coefficients  or  the  lattice  filter  reflection 

coefficients).  To  accurately  estimate  the  parameters,  the 

chosen  experimental  input  signal  must  be  sufficiently  rich 

in  frequency  content,  or  more  formally,  persistently 

exciting  so  that  it  excites  all  the  modes  of  the  system.  A 

sequence  is  said  to  be  persistently  exciting  of  order  n  if 

its  (n  x  n)  autocorrelation  matrix  is  nonsingular  [Ref. 

13:pp.  70-71].  By  using  the  techniques  presented  in  section 

D.3  of  this  chapter,  the  resulting  output  sequence  y(n)  can 

be  processed  to  yield  the  desired  reflection  coefficients. 

When  the  analysis  lattice  filter  is  implemented  with  these 

coefficients  embedded  in  its  structure,  the  lattice  becomes 

the  inverse  of  A  (z)  =  1/H(z).  This  procedure  is  depicted 

M 

in  Figure  2.10. 

The  application  of  linear  lattice  filters  to  decon¬ 
volution  was  simulated  using  the  FORTRAN  programs  LININV, 
ATOCOR ,  LEVIN,  AND  LATICE .  These  programs  are  provided  in 
Appendix  A.  The  simulation  was  conducted  for  a  second 
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order,  autoregressive,  LTI  system  defined  by  the  equation 
y ( n )  =  x ( n )  -  (0.6y(n-l)  +  0.08y(n-2))  .  As  previously 
mentioned,  modeling  the  system  was  the  first  step  in  the 
deconvolution  process.  In  order  to  identify  the  "unknown" 
parameters,  the  system  was  excited  by  a  zero  mean,  unity 
variance,  Gaussian  white  noise  sequence.  The  output 
sequence  y(n)  was  processed  by  subroutine  ATOCOR  to 
calculate  the  components  of  the  autocorrelation  matrix  R 

yy 

Then  subroutine  LEVIN  implemented  Levinson's  algorithm  to 
evaluate  both  the  prediction  error  filter  coefficients 
and  the  associated  reflection  coefficients  from  the  autocor¬ 
relation  matrix.  The  actual  output  of  subroutine  LEVIN  is 
the  lower  triangular  L  matrix  described  in  equation  (2.72); 
the  i-th  row  of  L  contains  the  i-th  order  prediction  error 
filter  coefficients  listed  in  reverse  order,  and  the  first 
column  contains  the  negative  of  the  reflection  coefficients. 
Once  the  reflection  coefficients  corresponding  to  1/H(z) 
were  calculated,  they  were  embedded  in  subroutine  LATICE 
which  implemented  the  analysis  lattice  filter  equations 
(2.76).  With  the  inverse  filter  of  H(z)  now  available 
(i.e.,  the  analysis  lattice  filter),  H(z)  was  driven  by  an 
"unknown"  signal  x(n).  The  output  of  H(z)  was  then  fed  into 
the  lattice  filter.  An  approximation  to  x(n)  was  recovered 
at  the  forward  error  output  signal  from  the  last  stage  of 


Simulations  were  run  both  with  and  without  the 


presence  of  measurement  noise.  Table  2.1  displays  the 
resulting  L  matrices  for  one  simulation.  The  case  of  no 
measurement  noise  is  shown  in  Figures  2.11  through  2.14. 
Figure  2.11  is  a  plot  of  the  input  signal  x(n),  Figure  2.12 
depicts  the  system  output  y(n),  and  Figure  2.13  shows  the 
lattice  filter  output  x(n).  As  can  be  seen,  the  results  of 
the  inverse  filtering  were  excellent — the  x  and  x  curves  are 
identical.  This  is  verified  by  Figure  2.14  which  shows  that 

A 

the  mean-square  error  between  x(n)  and  x(n)  is  nearly  zero. 
The  running  average  mean-square  error  was  calculated  using 
the  equation 


MSE(n) 


n 

(1/n)  V  (  x  (  i  ) 
i  =  l 


2 

x(i)  ) 


(2.97  ) 


The  simulation  was  then  repeated  for  a  nonzero 
measurement  noise.  Here,  the  added  measurement  noise  was  a 
zero  mean,  0.0025  variance,  Gaussian  white  noise  sequence. 
Using  the  same  input  as  for  the  previous  case,  the  outputs 
of  the  system  and  lattice  filter  are  shown  in  Figures  2.15 
and  2.16,  respectively.  Figure  2.17  is  a  plot  of  the  mean- 
square  error  between  x(n)  and  £(n).  The  results  in  Table  2.1 
show  that  even  with  the  presence  of  measurement  noise,  the 
system  parameters  were  still  accurately  identified.  The 
lattice  filter  recovered  the  basic  shape  of  x(n),  however 
it  could  not  remove  the  additive  measurement  noise 


Therefore,  the  output  of  the  lattice  filter  is  a  noisy 
or  "fuzzy"  version  of  x(n).  Additional  filtering  is  required 

A 

to  remove  the  noise  component  of  x(n).  The  results 
presented  here  are  representative  of  those  obtained  for 
simulations  involving  other  stable,  autoregressive,  LTI 


systems . 
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Figure  2.17  Mean-Square  Error  Between  x(n)  and  x(n) 
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Ill .  NONLINEAR  DECONVOLUTION 

A.  INTRODUCTION  TO  NONLINEAR  SYSTEMS 

While  the  previous  chapter  dealt  solely  with  linear 
systems  or  plants,  this  chapter  will  address  the  inverse 
filtering  problem  involving  nonlinear  systems.  The  chapter 
starts  with  a  brief  introduction  to  modeling  nonlinear 
systems.  Then  it  quickly  proceeds  to  extend  the  generalized 
linear  lattice  filter  results  of  section  II. D. 4  to  a 
nonlinear  analysis  lattice  filter.  The  nonlinear  lattice 
filter  is  discussed  in  detail.  To  conclude,  results  from 
numerous  simulations  involving  the  lattice  in  deconvolution 
applications  are  presented. 

System  linearity  is  defined  in  terms  of  the  principle  of 
superposition.  If  the  rule  by  which  the  system  transforms 
the  input  x(k)  into  the  output  y(k)  is  represented  by  the 
operator  T,  then  the  system  is  said  to  be  linear  if  and  only 
if 

T[ax  (k)  +  bx  (k)]  =  aT[x  (k)]  +  bT[x  (k)]  (3.1) 

12  1  2 

=  ay  (k)  +  by  (k) 

1  2 

for  arbitrary  constants  a  and  b.  The  system  is  nonlinear 
if  equation  (3.1)  is  not  satisfied.  Estimation  of  the 
parameters  of  a  nonlinear  system  is  a  complex  problem. 


One  approach  to  the  parameter  estimation  problem  for 
nonlinear  filters  is  the  Bayesian  approach.  This  approach 
leads  to  various  approximate  solutions  for  the  parameter 


estimates.  In  the  Bayesian  approach,  the  parameters  are 
considered  to  be  random  variables.  The  parameter  vector  h 
is  considered  a  random  vector  with  a  probability  density 
function  p(h).  Introducing  the  observation  vector  ydt)  and 
the  input  vector  x(  t )  ,  both  of  which  contain  data  up  to 
time  t,  the  a  posteriori  probability  density  function  for  h 
is  p  (  h  |  y_{  t )  ,  x  ( t )  )  .  One  possible  choice  for  the  estimate  of 

A 

the  h  vector  is  to  select  the  conditional  mean  h( t )  = 

E  [  h  |  yd  t )  j_x  (  t )  ]  .  Selecting  this  as  the  estimate  minimizes  the 

variance  of  the  parameter  estimation  error.  Another 

A 

possible  choice  is  to  select  the  h(t)  which  maximizes 
p(  h  |  y;(  t )  ,  x  ( t )  )  .  This  most  likely  value  is  known  as  the 
maximum  a  posteriori  (MAP)  estimate.  Finally,  the  Bayesian 
approach  also  leads  to  the  extended  Kalman  filter  for 

nonlinear  state  estimation  problems.  [Ref.  13:pp.  32-41] 

Another  popular  method  for  modeling  nonlinear  systems  is 
based  on  the  Volterra  series.  The  series  is  named  after  the 

mathematician  Vito  VOLTERRA.  The  first  person  to  apply  the 

series  to  nonlinear  systems  was  Norbert  WIENER.  If  the 
"black  box"  approach  is  taken  towards  a  nonlinear,  time- 
invariant  system,  the  relationship  between  the  input  x(t) 
and  the  output  y ( t )  can  be  represented  by  the  Volterra 
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where  n  =  1,2,...  and  h  (T  , . . . ,T  )  =  0  for  any  T  <  0, 

n  1  n  j 

j  =  l,2,...,n  .  The  functions  h  (T  ,...,T  )  are  called 

n  1  n 

Volterra  kernels  of  the  system.  This  equation  is  a 

functional  series.  That  is,  it  performs  an  operation  on  the 

function  x(t)  which  results  in  a  number  for  y(t).  If  the 

n-th  order  Volterra  operator,  H  ,  is  introduced  where 
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a 
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> 
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y  (t)  =  H  [ x ( t ) ] 
n  n 


rco  f-co 

■L 


(3.3) 


,  . . . ,T  )x(t-T  )...x(t-T  )  dT  . .  .dT  , 
n  1  n  1  n 


-oo  <-co 

then  the  series  can  be  expressed  in  operator  notation  as 


y ( t )  =  H  [x(t)]  +  H  [x(t)]  +  ...  +  H  [ x ( t )  ]  +...  (3.4) 

1  2  n 

oo  co 

=  y  H  [ x ( t ) ]  =  Y  y(t). 
n=l  n  n=l  n 
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Figure  3.1  is  a  graphic  representaion  of  this  equation. 


[Ref.  14:pp.  7-9]  The  H  operator  is  a  linear  operator, 

1 

H  is  a  quadratic  operator;  and,  in  general,  the  H  operator 
2  n 

involves  the  term  x(t)  to  the  n-th  power. 


B.  GENERALIZED  NONLINEAR  LATTICE  FILTER 

1 .  Introduction 

In  this  section,  it  will  be  demonstrated  how  the 
generalized  lattice  filter  introduced  in  section  II. D. 4  can 
be  applied  to  modeling  nonlinear  systems.  The  development 
will  start  by  looking  at  the  discrete  form  of  the  Volterra 
series.  Based  upon  this  series  an  alternate  tensor  notation 
representation  will  be  introduced.  The  results  of  section 
II. D. 4  will  then  be  extended  to  handle  a  two-dimensional 
field  of  data  which  will  result  in  the  generalized  nonlinear 
lattice  filter. 

2 .  Nonlinear  Lattice  Filter  Development 

The  discrete  form  of  the  Volterra  series  of  equation 
(3.2)  is  given  by 


y(k)  =  h  +  h  (n  )x(k-n  )  (3.5) 

0  11  1 

+  h  (n  ,n  )x(k-n  )x(k-n  )  +  ... 

2  12  1  2 


Using  LENK ’ s  tensor  notation,  an  equivalent  form  of  equation 


(3.5)  is  given  by 


■  "  wm  ■»«  w.  *  "WTr 


y ( k )  =  H  +  H  x  +  H  x  x  +  ...  ,  where  <3. 6a) 

Ai^j. 

A*  T 

x  '(  k  )  =  [  x  ( k  ),  x  (  k-1 ),...,  x  (  k- Xj ),...,  x  (  k-N  )  ]  (3.6b) 

for  A*,  =  0,1,..., N.  As  an  example,  if  N=1  and  equation 
(3.6a)  is  truncated  to  the  three  terms  shown  above,  then 
this  equation  can  be  rewritten  as 


y ( k )  =  H  +  H  x ( k )  +  H  x(k-l)  +  H  x(k)x(k) 

0  1  00 

+H  x(k)x(k-l)  +  H  x(k-l)x(k) 

01  10 

+  H  x( k-1 )x( k-1 )  .  (3.7) 

11 

Based  on  this  tensor  notation,  LENK  introduced  an  alternate 

representation  for  the  nonlinear  system.  Instead  of 

defining  the  components  of  the  vector  x  (k)  as  the  present 

and  past  values  of  the  signal  x(k)  as  in  equation  (3.6b), 

the  components  are  redefined  in  terms  of  x(k)  raised  to  the 

Aj  power  forA;  =  0,1,2,...,N.  That  is 

0  12  NT 

x'(k)  =  [x  (k),x  (k),x  ( k ),..., x  (k)J  (3.8) 

2  NT 

=  [  1  ,x(k)  ,x  ( k ),..., x  (k)] 


for  A*,  =  0,1,2,...,N.  Now  the  output  of  the  nonlinear  system 
is  defined  by 


f 

\ 
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for  K 

of  the 
term 


filter 

hva 


kernel;  it 
clarifly  the 
with  N=1  and 


=  0,1,2,.. .,P.  The  variable  P  gives  the  order 
and  N  defines  the  filter’s  finite  memory.  The 
plays  a  role  similar  to  that  of  the  Volterra 

N 

can  be  considered  a  (P+l)-order  tensor.  To 
meaning  of  this  notation,  the  following  example 
P=2  is  provided: 


y(k)  =  H.  1x  °(k)x  '(k-l)  . 

Aj  2 

=  H  +H  x(k)+H  x(k-l)+H  x(k)x(k-l)+H  x  (k) 
00  10  01  11  20 
2  2  2 
+  H  x  (k-l)+H  x  (k)x(k-l)+H  x(k)x  (k-l) 

02  21  12 


2  2 

+  H  x  (k)x  (k-l)  .  (3.10) 

22 

[Ref.  1 2 : pp .  39-48] 

The  above  nonlinear  model  for  y(k)  is  a  moving 
average  ( MA )  model:  the  output  is  defined  in  terms  of  the 

input  signal.  The  next  step  is  to  establish  an 
autoregressive  ( AR )  model  where  y(k)  is  estimated  in  terms 
of  its  past  values.  This  type  of  model  is  useful  when  the 
input  signal  is  not  readily  available.  An  autoregressive 
model  is  obtained  by  redefining  the  observation  vector  in 
terras  of  y(k)  vice  x(k).  Then  the  AR  model  is  given  by 


for 


y(k)  = 


=  1,2 


At  An 

y  (k-l)...y  (k-N)H  •  •  • -, 

A,  An 


.  ..,P.  If  the  system  is  driven  by 


(3.11) 


a  white 


noise  signal  u(k),  and  if  the  system’s  parameters  are 


approximated  by  H  .  ....  ,  then  the  estimation  error  is 

Ao  A** 

e(k)=y(k)-y(k)  =  [yA(k-l )  .  .  .yAt|k-N)H  ...  +  u(k)] 

A*  An 

-  [  yA‘(k-l  )  .  .  .y^k-NJH.  .  .  ]  (3.12) 

A|  ''N 

If  the  system  parameters  are  known  exactly,  then  the  estima¬ 
tion  error  and  the  white  noise  sequences  are  equal--that  is 
e(k)=u(k).  One  note  concerning  the  nonlinear  AR  model: 
unlike  the  linear  AR  model  whose  stability  is  easily 
determined  (i.e.,  the  system  is  stable  if  all  its  poles  lie 
within  the  unit  circle  in  the  complex  z-plane),  it  is 
difficult  to  judge  for  which  class  of  inputs  the  AR 
nonlinear  system’s  output  will  remain  bounded.  This  is 
because  the  order  of  the  nonlinearity  increases  with  time. 
[Ref.  1 2 : pp .  68-69] 

Using  this  alternate  tensor  form  of  the  AR  nonlinear 
system,  along  with  the  generalized  lattice  of  section 
II. D. 4,  the  nonlinear  lattice  filter  will  be  developed.  To 
simplify  the  discussion,  the  filter’s  finite  memory  will  be 
restricted  to  N  =  2.  Then  the  product  Y(k)  =  [ yA’( k-1  ) yA*( k-2  )  ] 
forAi,^j_  =0,1,..., P  forms  a  second  order  tensor,  or  a  two- 
dimensional  data  field  given  by: 


1 


y{k-2) 


p 

.  . .  y  (k-2) 


Y(k) 


P 


y(k-l) 

y(k-l ) y ( k-2  )  ... 

y ( k-1 ) y  (k-2) 

2 

2 

2  P 

y  (k-1) 

• 

y  (k-1 ) y ( k-2 ) . . . 

y  (k-l)y  (k-2) 

• 

• 

• 

P 

• 

• 

P 

• 

• 

P  P 

y  (k-1) 

y  (k-1 ) y (k-2 ) . .  . 

y  (k-l)y  (k-2) 

The  types  of  systems  that  this  nonlinear  lattice 
structure  can  model  exactly  are  of  the  form  shown  in  Figure 
3.2  .  The  nonlinear  combinations  block  forms  a  weighted  sum 
of  the  cross-products  and  powers  of  the  input  variables 
y(k-i),  for  i=l,2,...,N.  As  an  example,  with  N=2  and  P=2  the 
general  equation  for  the  system's  output  when  excited  by  the 
input  x(k)  is  given  by: 


y ( k )  =  x ( k )  -  {  H  +  H  y(k-l)  +  H  y(k-2) 

11  21  12 
2  2 

+  H  y(k-l )y(k-2 )  +  H  y  (k-1)  +  H  y  (k-l)y(k-2) 

22  31  32 

2  2 
+  H  y  (k-2)  +  H  y ( k- 1 ) y  (k-2) 

13  23 

2  2 

+  H  y  (k-l)y  (k-2)  }  (3.14) 

33 


It  is  also  assumed  that  the  system  is  time-invariant; 

otherwise  the  model  changes  with  time  k.  In  order  to  define 

2 

the  forward  and  backward  prediction  errors,  the  (P+1) 
elements  of  the  two-dimensional  data  field  must  by  converted 
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into  a  one-dimensional  sequence.  The  ordering  chosen  by 
LENK  is  listed  below: 

order  =  {( P, P ),( P- 1 , P ),< P , P-1 ),...,( 0 , P ),( P , 0 ) , 

(P-1,  P-1 ), (P-2,  P-1) , (P-1,  P-2), . . . , (0, P-1) , 
(P-1, 0) . (1,1), (0,1), (1,0), (0,0)}  ,  (3.15) 

where  the  components  of  the  data  field  are  identified  by  the 

indices  (ra,n)  for  m,n  =  0,1,..., P.  The  elements  of  the 

2 

ordered  set  are  numbered  consecutively  from  0  to  (P+1)  -1. 

The  notation  (m,n)-q  is  used  to  identify  the  q-th  element 

prior  to  the  element  (m,n)  as  referenced  to  the  ordered 

sequence.  This  notation  will  be  further  abbreviated  to 

mn-q.  Now  the  (q-l)-order,  normalized,  forward  error 

associated  with  predicting  the  value  of  the  element 
m  n 

y  (k-l)y  (k-2)  from  the  preceding  (q-1)  elements  is  given  by 

_q-i  q-1  A,  A, 

e  =  a  (m,n)y  (k-l)y  (k-2)  (3.16) 

mn  A)  Aj. 

v  <1'1 

for  At  .Aj.  =0 , 1 , 2  ,  .  .  .  ,  P .  Note  that  the  coefficients  a,  can 

A\ 

be  thought  of  as  the  components  of  a  second  order  tensor. 

q-1 

The  coefficient  a,  ■»  is  equal  to  zero  when  the  indices 

A\*x. 

(At  ,  A^)  do  n°t  correspond  to  the  (q-1)  elements  preceding 
(m,n)  in  the  ordered  sequence  (i.e.j  when  (  A » ,  )  >  (m,n) 

or  when  (Aj.Aj.)  <.  mn-q  ).  Also,  when  (A,»A^)  =  (m,n),  then 
q-1  q-1 

a  (m,n)=l/||e  it  .  (3.17) 

mn  mn 
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Similarly,  the  normalized  backward  prediction  error  in 
estimating  y(mn-q)  from  the  next  (q-1)  points  is  given  by 


_q-i  q-i  A,  Ax 

r  =  b  (mn-q)y  (k-l)y  (k-2) 
mn-q 


(3.18) 


q-1 


for  A,  ,Xx  =0,1,2, ... ,P.  The  coefficients  b.  ,  (mn-q)  equal 
zero  when  the  indices  (  A(,  A^)  do  not  correspond  to  the  (q-1) 
elements  following  (m,n)  (that  is,  when  (AiiAj.)  <  (mn-q)  or 
(  A,.  Al  )  >  (m.n)  ).  When  (Ai»Ai>  =  (“»n),  then 


q-1  q-1 

b  (mn-q)  =  1/  | I r  I  I . 

mn-q  mn-q 


(3.19) 


[Ref.  12 : pp.  145-146] 

Using  the  normalized  nonlinear  Levinson  algorithm, 
LENK  shows  that  the  nonlinear  prediction  errors  can  be 
updated  in  order  through  the  recursion  relation 


__q 

_q-l 
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e 

mn 

mn 

mn 

it 

_q 

q 

_q- 1 

r 

r 

mn-q 

mn-q 

(3.20) 


where 


mn 


mn 


-  K 


mn 


-K 


q 

1 


(3.21  ) 


and  the  unique  reflection  coefficients  are  given  by 


mn  q-1  q-1 

K  =  E{e  r  } 
q  mn  mn-q 


<3.22) 


[Ref.  12 :pp.  146-147] 


A  deficiency  remains  to  be  corrected:  the  goal  is 


to  estimate  y(k),  but  there  is  no  y(k)  term  in  the  two¬ 


dimensional  data  matrix  of  equation  (3.13).  This  problem  is 


solved  by  adding  another  channel  to  the  lattice  structure 


and  by  exploiting  the  orthogonality  of  the  backward  pre¬ 


diction  errors.  Since  the  backward  prediction  errors  leaving 


the  top  row  of  the  lattice  filter  (Figure  2.9)  are  uncor¬ 


related,  they  can  be  used  in  a  Fourier  series  to  estimate 


y(k).  These  backward  prediction  errors  can  be  formed  into  a 

2 

length  L  =  (P+1)  vector  defined  as 


X  0  12  _L-1  T 

[r  ]  =  [r  ,  r  ,r  , . . . ,r  ]  (3.23) 

(m,n)-/  00  00-1  00-2  00-L+1 


where  the  subscript  is  in  the  form  mn-q.  Now,  the  error  in 


estimating  y(k)  using  the  data  in  the  Y(k)  tensor  is 


calculated  from 


L  L-l  X 

e  =  y(k)  -  V  K,y  , 
k  A'  =0  A 


( 3.24  ) 


where’  the  Fourier  coefficients,  K  /  ,  are  given  by 
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[KJ  =  [E{y(k)r  },E{y(k)7  }  ,  •  •  •  ,  E  { y  (  k )  r  }  ]  .  (  3 . 25  ) 

A  00  00-1  00-L+1 

The  resulting  generalized  nonlinear  analysis  lattice 

structure  is  depicted  in  Figure  3.3.  [Ref.  12ipp.  147-150] 

The  nonlinear  lattice  structure  will  now  be  examined 

more  closely.  As  can  be  seen  in  Figure  3.3  ,  the  inputs  to 

the  analysis  filter  are  the  normalized  values  of  y(k)  and 
2 

the  (P+1)  +1  ordered  components  of  the  two-dimensional  data 

field  (equation  (3.13)).  The  ordering  of  these  inputs  is  in 

accordance  with  equation  (3.15).  At  a  given  time  k,  these 
2 

(P+1)  +2  inputs  are  evaluated  and  inserted  into  the  left 
side  of  the  lattice  structure.  The  lattice  calculations  are 
conducted  by  starting  at  the  top  row  and  moving  downward 
along  the  northwest-southeast  diagonal,  and  then  advancing 
to  the  top  of  the  next  diagonal  and  so  on.  When  all  the 
lattice  calculations  for  time  k  have  been  completed,  the 
process  is  repeated  for  time  k+1.  Just  as  in  the  case  of 
the  linear  lattice  filter,  the  PARCOR  coefficients  at  each 
lattice  section  act  to  decorrelate  the  two  input  signals. 
The  backward  error  signals  exit  at  the  top  of  the  lattice 
structure,  and  the  forward  error  signals  exit  at  the  right 
side.  The  output  of  interest  for  inverse  filtering  is  the 
forward  error  signal  from  the  top  row  of  the  filter  (i.e., 
the  row  corresponding  to  the  y(k)  input  signal).  This  is 
the  error  arising  from  the  estimation  of  y(k). 


O  1 


y  (k— l)y  (k— 2) 


3 .  The  Nonlinear  Lattice  Applied  to  Deconvolution 


¥ 


k 
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The  generalized  nonlinear  lattice  filter  can  be 
applied  to  the  inverse  filtering  problem  just  as  the  linear 
lattice  filter  was  in  section  II. D. 6.  The  first  problem  is 
one  of  system  identification  and  parameter  estimation.  Once 
the  model  parameters  have  been  estimated,  they  are  embedded 
in  the  nonlinear  filter  lattice  structure.  Then,  the 
nonlinear  lattice  represents  the  inverse  of  the  system’s 
transfer  function.  As  will  be  shown,  the  nonlinear  lattice 
filter  can  model  both  linear  and  nonlinear  systems--  the 
linear  system  is  just  a  special  case  of  the  nonlinear 
system.  This  inverse  filtering  algorithm  was  implemented  by 
the  FORTRAN  programs  NLMAIN ,  NLCLAT ,  SCHUR,  NORMS,  NLLAT , 
and  URAND .  These  programs  are  listed  in  Appendix  B.  Now, 
this  inverse  filtering  procedure  will  be  described  in  more 
detail . 

In  order  to  identify  the  model  parameters  (i.e.,  the 
PARCOR  coefficients),  the  system  was  excited  with  zero  mean, 
unit  variance  white  noise  sequences.  Both  Guassian  and 
uniform  noise  distributions  were  used  with  good  results. 
One  difficulty  encountered  in  generating  the  nonlinear  data 
was  ensuring  that  the  output  of  the  postulated  nonlinear 
system  remained  bounded  for  the  input  noise  signal.  For  the 
simulations,  the  filter  and  system  were  constrained  to  a 
finite  memory  of  at  most  two  delays  (N=2),  while  the 
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nonlinear  order  was  allowed  to  vary  from  P  =  0  to  P=4.  (For  N 


greater  than  two,  the  problem  of  ordering  the  elements  of 

the  N-dimensional  Y  tensor  increases  in  complexity.)  The 

system  output  sequence  y(k)  was  then  processed  by 

subroutines  NLCLAT  and  SCHUR,  which  determined  the 

corresponding  autocorrelation  matrix  and  the  partial 

correlation  coefficients,  respectively.  In  an  effort  to 

improve  the  accuracy  of  these  calculations ,  noise  sequences 

of  up  to  5,000  points  were  used.  Additionally,  the  PARCOR 

coefficients  were  averaged  over  as  many  as  50  realizations 

of  the  input  noise  random  process.  Through  trial  and  error, 

it  was  found  that  the  best  inverse  filtering  results  were 

obtained  when  the  resulting  reflection  coefficients  were 

truncated  to  two  decimal  places.  The  output  of  subroutine 

2  2 

SCHUR  is  a  ((P+1)  +1  x  (P+1)  +1)  upper  triangular  matrix  of 
reflection  coefficients.  This  matrix  is  simply  overlaid 
atop  the  upper  triangular  shaped  nonlinear  lattice  structure 
to  place  the  reflection  coefficients  at  the  correct  filter 
sections . 

With  the  reflection  coefficients  calculated  and 

embedded  in  the  filter,  we  were  able  to  use  the  lattice  in 

an  inverse  filtering  application.  To  recover  the  input 

signal  x ( k )  ,  the  system’s  output  signal  y(k)  was  processed 

by  subroutines  NORMS  and  NLLAT .  The  function  of  NORMS  was 

2 

the  (P+1)  +1  input  signals  into  the  lattice. 


a  a 


to  normalize 


Then  subroutine  NLLAT ,  using  the  previously  calculated 
reflection  coefficients,  implemented  the  lattice  structure 
and  carried  out  the  lattice  filter  calculations.  The  forward 
error  signal  out  of  the  top  row  of  the  lattice  yielded  the 
normalized  estimate  of  x(k).  Since  the  inputs  to  the  lat¬ 
tice  filter  are  all  normalized,  the  outputs  must  be 
denormalized .  Therefore,  since  the  input  into  the  top  row  of 
the  lattice  is  divided  by  the  norm  of  y(k),  the  forward 
error  signal  from  the  top  row  of  the  lattice  is  multiplied 
by  the  norm  of  y(k).  (Note  that  if  y(k)  is  a  zero  mean 
sequence,  then  this  norm  is  equivalent  to  its  standard 
deviation.)  Also,  it  was  found  that  to  achieve  a  good 
estimate,  it  was  necessary  to  further  scale  this  error 
signal  by  dividing  it  by  the  norm  of  the  noise  generated 
sequence  y ( k ) . 

In  —order  to  verify  the  accuracy  of  the  reflection 

coefficients,  the  noise  generated  output  signal,  y(k),  was 

passed  through  the  inverse  lattice  filter  to  see  how  well 

the  filter  whitened  this  sequence.  The  effectiveness  of  the 

whitening  filter  was  evaluated  by  examining  the  mean-square 

error  (MSE)  between  the  white  noise  input  signal,  x(k),  and 

A 

the  lattice  filter’s  output,  x(k).  The  running  average 
mean-square  error  was  calculated  using  the  equation 

VI  “  -  2 

( 1/k)  £  ( x ( i )  —  x ( i ) )  .  (3.26) 
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plot  essentially  provides  a  percentage  error  between  x(k) 
and  x(k).  Plots  of  the  MSE  are  included  in  the 
simulation  results.  Having  demonstrated  that  the  lattice 
filter  represented  a  good  inverse  of  the  system  H(z),  the 
system  was  then  driven  by  a  known  input  signal  x(k).  To 
recover  an  approximation  to  this  signal,  the  corresponding 
system  output  was  passed  through  the  lattice  filter,  and  the 
resulting  lattice  filter  output  was  denormalized  and 
rescaled  as  previously  discussed.  Simulation  results  for 
various  systems  are  presented  in  the  following  section. 

4 .  Inverse  Filtering  Simulation  Results 

In  this  section,  the  previous  modeling  and  inverse 
filtering  procedures  are  implemented  and  applied  to  various 
linear  and  nonlinear  systems.  Unless  otherwise  noted,  the 
reflection  coefficients  for  each  of  the  systems  were 
determined  by  using  twenty-five  realizations  (5,000  points 
each  )  of  the  zero  mean,  unit  variance  white  noise  random 
process  as  the  input  excitation  signal.  As  previously 
mentioned,  the  mean-square  error  between  this  white  noise 
input  and  the  lattice  filter's  output  was  plotted  to 
evaluate  the  lattice  filter’s  inverse  filtering  performance. 
After  the  system  was  modeled,  it  was  driven  by  the  signal 
x(k)  which  consisted  of  ramps,  pulses,  and  sinusoids  as 


shown  in  Figure  3.4. 


a.  System  I:  y{k)  =  x(k)  -  (0.6y(k-l)  +0.08y(k-2)) 

This  is  the  linear  system  first  introduced  in 

Section  II. D. 6.  An  input  Gaussian  white  noise  sequence  was 

used  to  model  the  system.  Figure  3.5  is  a  plot  of  the 

running  average  mean-square  error  between  the  input  noise 

and  output  error  signals.  In  Section  II. D. 6,  the  linear 

lattice  filter's  reflection  coefficients  were  found  to  be 

K  =  -0.555  and  K  =  -0.077.  It  should  be  noted  that  these 
1  2 

values  appear  in  the  first  row  of  the  nonlinear  lattice’s 
reflection  coefficient  matrix.  Here,  for  a  first  order, 
nonlinear  lattice  filter,  the  reflection  coefficients  are 
given  by  the  upper  triangular  matrix 

0  0  -.55  -.07  0 

0  0  0  0  -.43 

K  =  0  0  0  -.55  0 

0  0  0  0  0 

0  0  0  0  0 


The  system  output,  y(k),  corresponding  to  the  input  signal 
of  Figure  3.4  is  shown  in  Figure  3.6.  As  can  be  seen  by  the 
plots  of  x(k)  and  x(k)  in  Figure  3.7,  the  nonlinear  lattice 
filter  did  an  outstanding  job  of  recovering  the  "unknown" 
input  signal  x(k)  from  the  linear  system’s  output  y(k). 

b.  System  II:  y(k)  =  x(k)  -  ( 0 . 2y ( k-1 ) y ( k-2 ) ) 

This  system  involves  a  "cross-talk" 
nonlinearity,  that  is,  the  output  is  a  function  of  the 
product  of  two  different  signals.  For  this  system,  both 
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Figure  3.5  System  I  Mean-Square  Error 
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Gaussian  and  uniform  noise  distributions  were  used  in  the 
modeling  process  in  order  to  compare  the  two  techniques. 
Here,  both  methods  yielded  the  same  reflection  coefficients: 

0  0  0  0  -.21 
0  0  0  0  0 

K  =  0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

The  only  difference  between  the  two  techniques  was  the  value 
of  the  scaling  factor  (i.e.,  the  norm  of  the  noise  generated 
system  output  y(k)).  The  plot  of  y(k)  is  shown  in  Figure 
3.8.  The  running  average  mean-square  error  and  x(k)  plots 
for  the  Gaussian  noise  derived  model  are  depicted  in  Figures 
3.9  and  3.10,  respectively.  The  corresponding  plots  for  the 
uniform  noise  derived  model  are  given  in  Figures  3.11  and 
3.12.  As  can  be  seen  by  comparing  the  plots,  both  modeling 
variations  lead  to  nearly  identical  results.  Both 

techniques  yielded  excellent  approximations  to  the  input 
signal  x ( k ) . 

c.  System  III:  y(k)  =  x(k)  -  0 . 2y ( k- 1 ) y ( k- 1 ) 

This  system  involves  a  quadratic  nonlinearity. 
Therefore,  a  second  order  model  was  used.  The  system  output 
would  not  remain  bounded  for  a  Gaussian  noise  input,  so  the 
system  was  modeled  using  a  unit  variance  uniform  noise 
random  process.  The  resulting  reflection  coefficients  are: 
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Figure  3.9  System  II  Mean-Square  Error 
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Figure  3.11  System  II  Mean-Square  Error 
(Uniform  Noise  Model) 
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Figure  3.13  provides  the  running  average  mean-square  error 
plot.  The  output  y(k)  of  the  nonlinear  system  is  shown  in 


Figure  3.14.  Figure  3.15  depicts  the  comparison  between  the 

.  A  .  . 


system  input  x(k)  and  the  lattice  filter  output  x(k).  The 
x(k)  curve  has  the  same  shape  as  x(k),  however,  it  is 
slightly  offset.  In  an  attempt  to  improve  the  approximation 
process,  the  number  of  realizations  of  the  noise  random 
process  used  to  model  the  system  was  doubled  from  twenty- 
five  to  fifty  and  the  simulation  was  repeated.  Although 
several  reflection  coefficient  values  changed  by  .01,  there 
was  no  perceptible  improvement  in  the  estimation  of  x(k). 

d.  System  IV:  y(k)=x(k)  -  ( 0 . 5+0 . 6y ( k- 1 ) + . 08y ( k-2 ) ) 
This  example  consists  of  the  linear  System  I 
modified  by  the  addidtion  of  a  constant  bias  term.  Gaussian 
noise  was  used  to  model  the  system.  The  reflection 
coefficients  determined  for  the  first  order  model  are: 


.  A  A-W. 

*  M  .  .H  ■ 


ER 


mm 


•4 


55  -.07 
23  -.40  -. 


0  - , 

0 

0 


This  is  basically  the  same  matrix  as  obtained  for  System  I, 


but  modified  by  the  addition  of  several  terms  which  attempt 


to  model  the  constant  offset.  (Note  the  reappearance  of  the 


values  of  -.55  and  -.07  in  the  top  row  of  the  lattice 


matrix.)  The  running  average  mean-square  error  is  shown  in 


Figure  3.16.  The  system  output  y(k)  is  plotted  in  Figure 


3.17.  Figure  3.18  is  a  plot  of  the  lattice  output  x(k) 


and  the  system  input  x(k).  The  x(k)  curve  is  an  excellent 


replica  of  the  original  input  signal,  but  it  is  offset  by  a 


constant  value.  Although  the  small  mean-square  error  evident 


in  Figure  3.16  indicates  that  the  lattice  is  a  good  inverse 


filter,  the  lattice  was  unable  to  completely  remove  the  bias 


term.  Increasing  the  order  of  the  model  (i.e.,  overmodeling 


the  system)  did  not  improve  the  results. 


e.  System  V:  y(k)=x(k)  -  ( 5 . 0+0 . 6y ( k- 1 ) + . 08y ( k-2 ) ) 


In  order  to  further  investigate  the  constant 


offset  nonlinearity,  the  example  of  System  IV  was  repeated 


using  a  bias  of  5.0  instead  of  0.5.  Again,  Guassian  white 


noise  was  used  in  determining  the  model’s  parameters.  The 


reflection  coefficients  of  the  first  order  nonlinear  lattice 


are  given  by: 


OU3SW 

Figure  3.16  System  IV  Mean-Square  Error 
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K 


h* 


0  -.92  -.55  -.07  0 

0  •  0  -.92  -.86  -.74 

0  0  0  .78  -.75 

0  0  0  0  -.90 

0  0  0  0  0 


Again,  the  reflection  coefficient  values  of  -.55  and  -.07 
occur  in  the  top  row  of  the  matrix  as  these  elements 
apparently  model  the  linear  portion  of  the  sytem.  The 
running  average  mean-square  error  is  displayed  in  Figure 
3.19,  and  the  output  of  the  nonlinear  system  is  shown  in 
Figure  3.20.  Figure  3.21  provides  the  comparison  between  the 
system  input  and  the  output  of  the  inverse  filter.  As  with 
System  IV,  the  x(k)  curve  produced  by  the  deconvolution 
process  has  the  same  shape  as  x(k),  but  is  offset  from  the 
desired  curve  by  a  small  constant  value.  Here,  the  system 
bias  is  5.0,  and  the  x(k)  and  x(k)  curves  differ  by  about 
0.5.  This  is  a  relative  improvement  over  the  results 
obtained  for  System  IV.  System  IV  had  a  constant  bias  of 
only  0.5  but  the  x{k)  curve  was  offset  from  the  x(k)  curve 
by  approximately  0.4. 
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IV.  CONCLUSION 


In  this  thesis,  several  linear  least-squares 
deconvolution,  or  inverse  filtering,  techniques  were 
reviewed.  Particular  emphasis  was  placed  on  the  lattice 
filter.  It  was  shown  that  when  a  system  is  defined  by  an 
autoregressive,  linear  time-invariant  model,  this  model  can 
be  transformed  into  an  equivalent  lattice  filter 
representation.  Furthermore,  the  resulting  analysis  lattice 
filter  acts  as  a  whitening  filter,  and  is  the  inverse  of  the 
system’s  autoregressive  transfer  function.  An  example 
demonstrated  the  effectiveness  of  the  lattice  filter  in  a 
linear  deconvolution  application.  Unfortunately,  the  linear 
lattice  filter  is  unable  to  model  nonlinear  systems  and, 
therefore,  is  not  a  viable  nonlinear  inverse  filtering 
technique . 

The  discussion  of  the  linear  lattice  filter  led  to  the 
development  of  the  generalized  lattice  filter,  which  in  turn 
led  to  the  derivation  of  the  nonlinear  lattice  filter. 
The  goal  of  this  thesis  was  to  implement  the  nonlinear 
lattice  filter,  arid  then  apply  it  to  the  linear  and 
nonlinear  deconvolution  problem.  This  was  accomplished  with 
generally  very  good  results.  It  was  shown  that  the 
nonlinear  lattice  filter  was  suitable  for  modeling 


I 


discrete  autoregressive  systems  where  the  output  y(n) 

consisted  of  a  weighted  sum  of  the  cross-products  of 
P  q 

the  terms  y  (n-1)  and  y  (n-2),  where  the  powers  p  and  q  can 
take  on  the  values  {0,1, 2, 3, 4}.  Examples  were  presented 
demonstrating  the  ability  of  the  nonlinear  lattice  filter  to 
effectively  act  as  an  inverse  filter  for  both  linear  and 
nonlinear  autoregressive  systems. 
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APPENDIX  A 

LINEAR  LATTICE  FILTER  FORTRAN  PROGRAMS 


***************  ************************* ************************** 


LT  SCOT  L.  JOHNSON 
LINEAR  INVERSE 
06  APRIL  1986 


****************************************************************** 


THE 
SYSTEM 


PURPOSE  OF  THIS  PROGRAM  IS  TO  MODEL  A  LINEAR  AUTOREGRESSIVE 

_ TEM  BY  THE  EQUIVALENT  ANALYSIS  LATTICE  FILTER. 

THE  SYSTEM'S  TRANSFER  FUNCTION  IS  GIVEN  BY  H(Z)  =  1/A(Z) 

WHERE  A(Z)  =  (Z**2  +  A1*Z  +  A2)/Z**2. 

THE  LATTICE  FILTER  REFLECTION  COEFFICIENTS  ARE  DETERMINED 
BY  DRIVING  H(Z]  WITH  ZERO  MEAN,  UNITY  VARIANCE  GAUSSIAN  WHITE 
NOISE.  THE  OUTPUT  IS  PROCESSED  BY  SUBROUTINE  ATOCOR  WHICH  COMPUTES 
THE  COMPONENTS  OF  THE  SAMPLED  AUTOCORRELATION  MATRIX  R.  SUBROUTINE 
LEVIN  IMPLEMENTS  THE  LEVINSON  ALGORITHM  AND  EVALUATES  THE  LATTICE 
FILTER  COEFFICIENTS  FROM  THE  AUTOCORRELATION  MATRIX.  FINALLY, 
SUBROUTINE  LATICE  IMPLEMENTS  THE  ANALYSIS  STRUCTURE  OF  THE 
LATTICE  FILTER  WHICH  IS  EQUIVALENT  TO  THE  INVERSE  FILTER  H(Z). 

NOW  WHEN  H(Z)  IS  DRIVEN  BY  AN  UNKNOWN  SIGNAL  xtNj,  THE 
SYSTEM  OUTPUT  Y(N)  CAN  BE  FED  INTO  THE  PREVIOUSLY  DETERMINED 
LATTICE  FILTER:  SINCE  THE  LATTICE  FILTER  IS  THE  INVERSE  OF  H(Z) 

THE  LATTICE  FILTER  OUTPUT  SHOULD  BE  A  GOOD  ESTIMATE  OF  X(N). 

***VARIABLE  DEFINITIONS**** 

Al.  A2  =  COEFFICIENTS  OF  THE  PREDICTION  ERROR  POLYNOMIAL  A(Z) 
DSGED.DSEEDI  =  SEED  VALUES  USED  BY  THE  IMSL  WHITE  GAUSSIAN  NOISE 
GENERATOR  FUNCTION  GGNQG 

E  =  VECTOR  OF  MEAN-SQUARED  PREDICTION  ERRORS  E(0) ,E(1) . E(ORDER) 
ESUM  =  VECTOR  USED  IN  DETERMINING  AVERAGE  E 
GAMMA  =  VECTOR  OF  LATTICE  REFLECTION  COEFFICIENTS.  GAMMA(I)  IS 
THE  REFLECTION  COEFFICIENT  OF  THE  IYTH  LATTICE  STAGE. 
GAMSUM  =  VECTOR  USED  IN  DETERMING  AVERAGE  GAMMA 
GRAFP  =  REAL  VALUE  WHICH  DEFINES  LENGTH  OF  X-AXIS  FOR  PLOTTING 
L  =  LOWER-TRIANGULAR  MATRIX  WHOSE  ROWS  ARE  THE  REVERSE  OF  ALL  THE 
PREDICTION  ERROR  FILTERS  FROM  ORDER  ZERO  TO  THE  HIGHEST  ORDER 
LSUM  =  MATRIX  USED  IN  DETERMINING  THE  AVERAGE  L  MATRIX 
MSE  =  MEAN-SQUARE  ERROR  BETWEEN  X(N)  AND  XHAT(N) 

MSEMAX  =  MAXIMUM  VALUE  OF  MSE 

MSESUM  =  RUNNING  SUM  USED  IN  CALCULATING  MSE 

NINDEX  =  ARRAY  OF  REAL  NUMBERS  USED  IN  DISSPLA  PLOTTING  ROUTINES 

NOISE  =  ARRAY  OF  MEASUREMENT  NOISE  ADDED  TO  OUTPUT  OF  H(Z) 

NUMPTS  =  NUMBER  OF  POINTS  IN  INPUT  NOISE  SEQUENCE 
ORDER  =  ORDER  OF  THE  LATTICE  FILTER 
ORDERP  =  ORDER  +  1 

PLTPTS  =  NUMBER  OF  POINTS  USED  IN  PLOTTING  ROUTINES 
R  =  VECTOR  OF  AUTOCORRELATION  LAGS  R(0) ,R( 1)A. .  .  R(ORDER) 

RMAX  =  MAXIMUM  MAGNITUDE  OF  ELEMENTS  OF  Rl  TG  BE  USED  IN  PLOTTING 
STDDEV  =  STANDARD  DEVIATION  OF  THE  MEASUREMENT  NOISE 
TRIAL  =  NUMBER  OF  REALIZATIONS  OF  THE  WHITE  GAUSSIAN  NOISE  RANDOM 
PROCESS  USED  IN  MODELING  THE  SYSTEM  H(Z) 

X  =  INPUT  SEQUENCE  INTO  THE  SYSTEM  H(Z) 

XMAX,  XMIN  =  RANGE  OF  X  VALUES  USED  IN  DISSPLA  PLOTTING  ROUTINES 
XHAT  =  ESTIMATE  OF  X-  OUTPUT  OF  THE  ANALYSIS  LATTICE  FILTER  A(Z) 

Y  =  OUTPUT  SEQUENCE  $F  H(Z) 

YMAX ,  YMIN  =  RANGE  OF  Y  VALUES  USED  IN  DISSPLA  PLOTTING  ROUTINES 

***VARIA8LE  DECLARATIONS**** 

INTEGER  I, N,K, NUMPTS, ORDER, ORDERP, PLTPTS, TRIAL 
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-  ■  A. 
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10 


11 

c 


12 

C 

C 


18 

C 


REALM  X( 5000) ,Y( 5000) ,A1 ,A2,XMIN ,XMAX , YMAX.YMIN ,NINDEX(  5000) 
REALM  GAMMAfl6),GAMS0M(ld).ft(10)  L( 10  10) . £( 10) ! RMAX .GRAFP 
REALM  XHATf5O00j,LSUMno.i6);ES0fl(10)  N01§E(5006),STODEV 


PEAL* 4  MSE(5OO0)lMSESU&,M$EPl!.t 
REAL*8  DSEED.DSEfOl 


DEFINE  THE  LENGTH  OF  THE  TIME  SEQUENCES 
NUMPTS  =  5000 
PLTPTS  =  2000 
GRAFP  =  FL0AT(PLTPTS+1) 


DEFINE  THE  ORDER  OF  THE  LATTICE  FILTER  AND  THE  A(Z)  COEFFICIENTS 
ORDER=2 
0RDERP=0RDER+1 
A1  =  0.6 
A2  =  0. 08 

INITIALIZE  VARIABLES 
DO  6  1=1 ,0RDERP 
GAMMA( ’ ' 

GAMSUt*  , 

ESUM( I , 

Rin  =  o. 

DO  4  K=1 .ORDERP 
LSUM( I ,K)  =  0. 

CONTINUE 
CONTINUE 


JtKr 


XMAX  =  0. 

YMAX  =  0. 

MSESUM  =  0. 

MSEMAX  =  0. 

DEFINE  THE  SEEDS  FOR  THE  NOISE  SIGNALS,  THE  NUMBER  OF  NOISE 
REALIZATIONS  USED  IN  MODELING  H(Z),  ANO  THE  STANDARD  DEVIATION  OF 
THE  MEASUREMENT  NOISE 
DSEED  =  1243073. 5D0 
DSEED1  =  724389. 4D0 
TRIAL  =  25 
STDDEV  =  0.  05 


DO  25  1=1 .TRIAL 

DSEEti  =  DSEED/DFLOAT(I ) 

DSEED1  =  DSEED1/DFL0AT(I) 

DO  10  K=l, NUMPTS 

X(K)  =  GGNQF(DSEED) 

CONTINUE 

DO  11  K=l. NUMPTS 

NOIS£(K)  =  STDDEV  *  GGNQF(DSEEDl) 
CONTINUE 


ffl  :  5(iJ  - 

6  12  K  =  3 ,NL 

\  /  /  iy  \  ’w  / 


noise  m  ,  x 
A1*Y(1)  +  NO I S  E ( 2 ) 

.  .  _  .NUMPTS 

Y(K)  =  X(K)  -  (A1*Y(K-1)  +  A2*Y(K-2))  + 
CONTINUE 


NOISE(K) 


DETERMINE  THE  AUTOCORRELATION  VECTOR  OF  THE  SEQUENCE  Y(N) 
CALL  ATOCOR( NUMPTS, Y .ORDER, R.RMAX) 

DETERMINE  THE  L  MATRIX  AND  E  VICTOR 
CALL  LEVINToRDER.R.L.E) 

DO  18  K=l. ORDER 

GAMMMK)=-1*L(K+1,1) 

GAMSUM(K)  =  GAMSUM(K)  +  GAMMA(K) 

CONTINUE 


DO  22  N=l, ORDERP 

ESUM(N)  =  ESUM(N) 

DO  20  K=1 .ORDERP 

LSU^(N.K)  =  LSUM(N.K) 


E(N) 


L(N,K) 
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CONTINUE 

CONTINUE 

CONTINUE 


CALCULATE  THE  AVERAGE  VALUES  OF  E,  L,  AND  GAMMA 
DO  28  N=1  .ORDERP 

EfN)  =  ESUM(N)/FLOAT( TRIAL) 

DO  26  K=1 .ORDERP 

L(N,K)  =  L$UM(N,K)/FLOAT(TRIAL) 

CONTINUE 
CONTINUE 
DO  30  K=l, ORDER 

GAMMX(K)  =  GAMSUM(K)/FLOAT( TRIAL) 

CONTINUE 


WITH  INPUT  Y  AND  LATTICE  PARAMETERS  GAMMA.  DETERMINE  THE 
OUTPUT  OF  THE  ANALYSIS  MODEL  LATTICE  FILTER 
CALL  LATICE(NUMPTS, ORDER, GAMMA, Y.XHAT) 


PRINT  RESULTS 


PRINT  THE  R.E,  AND  GAMMA  VECTORS 
WRITE(8.35) 

FORMAT fTi, ' 1 ' ,T4, 1 N* , T10 ,  R( N) 1 ,T20, 1 E(N) 1 ,T30, 'GAMMA(N)  ) 
NULL-U 


NULL=D 

WRITE(8 , 40)  NULL. R(l).Eif  1) 

FORMAT(Tl  'O'  J4,f2,Tld,3F10.4) 

DO  45  K=2,OR£oR 
N=K-1 

WRITE(8,40)  N,R(K) ,E(K) ,GAMMA(N) 
CONTINUE  ' 


PRINT  THE  L(I,J)  MATRI 
WRITE(8,48) 

FORMATCTi'I1  'L(I,J)= 
DO  55  I=LORt)ERP 


MATRIX 


FO®Ttfff:?00Uo^FIiOJ)4Sx15)MDERP) 

CONTINUE  1 


DO  THE  DECONVOLUTION  PROBLEM.  FIRST  DEFINE  THE  DETERMINISTIC 
INPUT  X(N),  THEN  DETERMINE  THE  SYSTEM  OUTPUT. 

DO  85  K=1  PLTPTS 

X(K)  =  1.0  *  SIN(0. 0 1 2 6 * F LOAT ( K ) ) 

CONTINUE 


(1)  =  XI 

)  +  NOISEf 

L 2 5  =  XC2 

)  -  A1*Y(1 

10  90  K=3.l 

^LTPTS 

NO I S  E( 2 ) 


Y(K)  =  X(K)  -  (A1*Y(K-1)  +  A2*Y(K-2))  +  NOISE(K) 
CONTINUE  ' 


INPUT  Y(N)  INTO  THE  LATTICE  TO  RECOVER  THE  ESTIMATE  OF  X(N) 
—  '  LATICEr'  ~ . . 


PLTPTS, ORDER, GAMMA, Y.XHAT) 


CALCULATE  THE  MEAN-SQUARE  ERROR  (MSE)  BETWEEN  X(N)  AND  XHAT(N) 
DO  93  K=1 APLTPTS 

MSE$UM  =  [ X( K)-XHAT(K))*(X( K)-XHAT( K) )  +  MSESUM 
MSEQK)  =  SQRT(MSESUM/FLOAT(K)) 

IF(MSE(K). GT. MSEMAX)  MSEMAX  =  MSE(K) 

CONTINUE 


DEFINE  PLOTTING  PARAMETERS 
DO  95  K=  KPLTPTS 


XMAX  =  ABS(XfK)) 

AX)  XMAX  =  ABm  XHAT(K)) 
YMAX  =  ABS( Y( K)) 


.‘  V.  -ft,: ^ Vs o ;  rj»  >  *.»>.* 


& 

95 

CONTINUE 

ft*. 

XMIN  =  -1.  0 

*  XMAX 

h 

r 

YMIN  =  -1.0 

*  YMAX 

k « 


C 

C 

c 

c 


c 

c 


c 

c 


****  DISSPU  PLOTTING  ROUTINES  **** 

PLOT  THE  SYSTEM  INPUT  AND  THE  LATTICE  OUTPUT 
CALL  TEK618 
CALL  COMPRS 
CALL  RESET ( ' ALL 1 ] 

CALL  PAGECo.  0,6. b) 

CALL  XINTAX 
CALL  AREA2DC6.  75,5.  0] 

CALL  XNAMEPTIME  $AMPLES$'  ,100) 

CALL  YNAMEC ' X(N)$r , 100) 

CALL  CROSS 

CALL  GRAFCO. ,500.  ,GRAFP. XMIN.  SCALE1 ,XMAX) 

CALL  CURVEfNiNDEX  X.PLT^TS.O) 

CALL  ENDPL(O) 

PLOT  THE  SYSTEM  OUTPUT  Y(N) 

CALL  NOBRDR 
CALL  AREA2D(6. 75,5.0] 

CALL  XNAME(‘TIME  SAMPLES$' ,100) 

CALL  YNAMEC 'Y(N)$  ,100] 

CALL  HEADINT ' YCN)  =  X(N)  -  (0. 6*Y(N-l)+0. 08*YfN-2))S'  100,1. 5,1] 
CALL  HEADINf 1  V(N]=X(N)-(0.  6*Y(N-1  )+0.  08*Yf  N-2) )+NOISi:$ ' ,100,1.  5,2) 
CALL  HEADINC 1  NOISE  =  N(0,0.  0025)$\  100, 1.  5,2) 

CALL  CROSS  ,  .  ’ 

CALL  GRAFfO. ,500.  .GRAFP.YMIN, ' SCALE1 ,YMAX) 

CALL  CURVEfNiNDEX  Y.PLTpTS.O) 

CALL  ENDPL(O) 

PLOT  LATTICE  FILTER  OUTPUT,  XHAT(N) 

CALL  AREA2D(6. 75,5.0] 

CALL  XNAMEPTIME  SAMPLESS'.IOO) 

CALL  YNAMEC 'XHAT(N)$' ,100) 

CALL  CROSS 

CALL  GRAFfO. ,500.  .GRAFP.XMIN, 1  SCALE' ,XMAX) 

CALL  CURVEfNiNDEX  XHAT.PLTPTS.O) 

CALL  ENDPL(O)  -  — 

PLOT  THE  MEAN-SQUARE  ERROR  BETWEEN  X(N)  AND  XHAT(N) 

CALL  AREA2D  6.75,5.0] 

CALL  XNAMEpTIME  SAMPLES^ ' ,  100) 

CALL  YNAMEC 1 MSE( N)$ ' ,100) 

CALL  CROSS 

C  CALL  GRAFfO.  ,500.  ,GRAFP  ,0. ,' SCALE ' ,MSEMAX) 

CALL  GRAFfO.  ,5001  ,GRAFf>,o!p SCALE'  ,0.08) 

CALL  CURVEfNiNDEX  MSE,P£.TPtS,0) 

CALL  ENDPL(O) 

CALL  DONEPL 
C 

STOP 

END 

C 

C 

£★**********★********★*********************★★****★****★*****★*#********* 

C 
C 
C 
C 

c 
c 
c 
c 


THE  SAMPLED 


SUBROUTINE  ATOCOR 

GIVEN  A  TIME  SEQUENCE  YfN],  THIS  PROGRAM  CALCULATES 
AUTOCORRELATION  MATRIX  TERMS  R(0) ,R(  1) ,.  .  .  , R( ORDER ) . 

WRITTEN  06  APRIL  1986 

£  ★  *****  ★  *  *  *  ★  ★★  *  *  *  *  *  *  ★  *  *  *  *  *  rt  *  *  *  ★  rt  *  *  *  ★  *  *  *  *  *  *  rt  *  *  *  *  *  *  *  *  *  W  rt  *  rt  *  i*t  *  ★  ***********  * 

C 


.  SUB ROUT INE^ ATOCORj NUMPTS , Y , ORDER , R , RMAX ) 


VARIABLE  DEFINITE 
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Y( 5000)=  INPUT  SEQUENCE 

Rf 10)  =  AUTOCORRELATION  MATRIX 

NOMPTS  =  NUMBER  OF  POINTS  IN  THE  SEQUENCE  Y(N) 

ORDER  =  MAXIMUM  LAG  FOR  WHICH  R  IS  EVALUATED 

OROERP  =  ORDER+1: THE  LENGTH  OF  THE  R  VECTOR 

SUM=  RUNNING  TOTAL  OF  PRODUCTS  FOR  A  GIVEN  LAG 

RMAX=  MAXIMUM  VALUE  OF  R;  USED  FOR  DISSPLA  PLOTING 

INTEGER  NUMPTS, I ,J ,K,N,ORDER,ORDERP 
REAL  Y(5000),R(l6),SClM  RMAX 
RMAX=0. 

ORDERP=ORDER+l 
DO  75  K=1 .ORDERP 
SUM=0. 

N=NUMPTS-K+1 
DO  70  J=1,N 

SUM=SUM+(Y(K+J-1)*Y( J)) 

CONTINUE 

R(K)=SUM/NUMPTS 

IF(ABS(R(K)). GT.  RMAX)  RMAX=ABS(R(K)) 

CONTINUE  v  v 

RETURN 


SUBROUTINE  LEVIN 

THIS  SUBROUTINE  IMPLEMENTS  LEVINSON'S  ALGORITHM.  IT  GENERATES  ALL 
THE  PREDICTION  ERROR  FILTERS  UP  TO  A  GIVEN  ORDER,  FROM  THE 
AUTOCORRELATION  LAGS. 

BASED  ON  A  PROGRAM  WRITTEN  BY  S.J.  ORFANIDIS  (REF.  1:  P.  333) 
MODIFIED  06  APRIL  1986 

*********************************************************************** 

SUBROUTINE  LEVI N( ORDER, R, L.E) 

***VARIABLE  DEFINITIONS’**^* 

ORDER  =  ORDER  OF  LATTICE  FILTER  '  ~ 

R  =  VECTOR  OF  AUTOCORRELATION  LAGS  Rf 0) ,R( 1),. . . ,R(ORDER) 

L  =  UNIT  LOWER-TRIANGULAR  MATRIX.  ITS  1-TH  ftOW  AoLDS  THE  I-TH 
PREDICTION-ERROR  FILTER  IN  REVERSE  ORDER.  ITS  FIRST  COLUMN 
HOLDS  THE  NEGATIVES  OF  THE  REFLECTION  COEFFICIENTS, 

GAMMA(I)  =  -L( 1+1,1)  FOR  1=1,2.. .. , ORDER 
E  =  VECTOR  OF  PREDICTION  ERRORS  Efo),E(l), _ E(ORDER) 


,E(ORDER 


E  =  VECTOR  OF  PREDICTION  ERRORS  E( 0)  , E( 1 ) ,. . . ,E(ORDER) 

THE  MATRIX  L  AND  THE  DIAGONAL  MATRIX  D  FOftMED  BY  THE  E'S  DEFINE 
A  UL  CHOLESKY  FACTORIZATION  OF  THE  INVERSE  OF  THE  AUTOCORRELATION 
MATRIX:  R  INVERSE  =  L  TRANSPOSE  *  D  INVERSE  *  L 


***VARIABLE  DECLARATIONS**** 

REAL  R( 10) .E( 10) . LC 10, 10) ,GAPaGAMMA 
INTEGER  r|ptUS;iMlNU$,j;flRDEft,ORDERP 
ORDERP=OR0ER+l 

SET  THE  UPPER  TRIANGLE  OF  THE  L  MATRIX  TO  ZERO 
DO  60  1=1, ORDER 
IPLUS=I+1 
IMINUS=I-1 

DO  60  J=IPLUS,ORDERP 
L(i,j)=d. 

CONTINUE 

Lf  1,11=1. 

LC2, 2)=1.0 

Lf  2  15=-R(2)/R(1) 

E  1  =R(1) 


§ 

§ 


<s 


3 


<Jd 

fjy 


GAP=0. 

IMINUS=I-1 
DO  63  K=1 , IMINUS 

GAP=GAP+R(K+1)*L(I-1,K) 

CONTINUE 

GAMMA=GAP/E( 1-1) 

L?I,1)=-1*GAMMA 
DO  64  K=2, IMINUS 

L(I  K)=L(I-1,K-1)-GAMMA*L(I-1,I-K) 

CONTINUE 


CONTINUE 

RETURN 

END 


1)*( 1-GAMMA**2) 


C 

C  SUBROUTINE  LATICE 
C 

C  THIS  PROGRAM  IMPLEMENTS  A  SINGLE  CHANNEL  LATTICE  STRUCTURE. 

C  WHEN  GIVEN  THE  LATTICE  COEFFICIENTS  AND  THE  INPUT  SEQUENCE, 

C  THE  PROGRAM  DETERMINES  THE  OUTPUT  SEQUENCE. 

C  WRITTEN  06  APRIL  1986 

C 

£*********************************************************************** 

C 

SUBROUTINE  LATICEf NUMPTS. ORDER, GAMMA, Y, OUTPUT) 

C  ***VARIABLE  DEFINITIONS**^ 

C  NUMPTS=  NUMBER  OF  POINTS  IN  THE  SEQUENCES;  MAX  IS  5000 
C  ORDER=  ORDER  OF  THE  LATTICE;  MAX  IS  9 
C  GAMMA( ORDER )=  LATTICE  COEFFlCENT  ARRAY 

C  F  =  FORWARD  ERROR 

C  B  =  BACKWARD  ERROR 

C  DELAY?ORDER)=  ARRAY  OF  DELAYED  BACKWARD  ERROR  SIGNALS 

C  TEMPf ORDER)  =  ARRAY  WHICH  TEMPORARILY  HOLDS  THE  BACKWARD  ERROR 
C  Y(NUMPTS)=lNPUT  DATA  ARRAY 

C  OUTPUT( NUMPTS )=ARRAY  OF  LATTICE  OUTPUT  DATA 

C  ***VARIABLE  DECLARATIONS**** 

INTEGER  NUMPTS .ORDER , I ,  K.M 

I REAL^GAMMA^ 10^ F,B,DfLAY( 10) ,TEMp(IO) ,Y(5000) ,OUTpuT( 5000) 

DO  80  1=1  .ORDER'* 

DELAY(I)=0. 

TEMP( I )=0. 

80  CONTINUE 

C 

C  DO  TIME  ITERATION 

DO  88  K=l, NUMPTS 
F=Y{K) 

B=Y( K) 

C  FOR  EACH  TIME  INSTANT, RECURSIVELY  INCREASE  THE  LATTICE  ORDER 

DO  85  M=1 .ORDER 
TEM^(M)=B 

B=  D  E  LA  Y  ( M )  - ( GAMMA ( M ) *  F ) 

F=F-(GAMMA(M)*DELAY(M) ) 

DELAY(M)=TEMP(M) 

85  CONTINUE 

OUTPUT(K)=F 
88  CONTINUE 
RETURN 
END 
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APPENDIX  B 

NONLINEAR  LATTICE  FILTER  FORTRAN  PROGRAMS 


*********************************************************************** 

* 

PROGRAM  NLMAIN 

* 

THIS  PROGRAM  DEFINES  THE  SYSTEM  PARAMETERS  AND  THE  SYSTEM'S 
INPUTS  AND  OUTPUTS.  IT  CALLS  SUBROUTINES  TO  DETERMINE  THE 
THE  CORRESPONDING  NONLINEAR  ANALYSIS  LATTICE  MODEL. 

* 

WRITTEN  30  APRIL  1986 

* 

*********************************************************************** 

THIS  PROGRAM  INPUTS  A  WHITE  NOISE  SEQUENCE  X(K)  INTO  A  AUTOREGRES¬ 
SIVE.  NONLINEAR  SYSTEM  H(Z)  .  THE  SYSTEMYS  OUTPUT,  Y(K),IS 
PROCESSED  TO  DETERMINE  THE  AUTOCORRELATION  MATRIX  AND  THE  NONLINEA 
LATTICE'S  REFLECTION  COEFFICIENT  MATRIX.  SINCE  THE  OUTPUT  OF  THE 
LATTICE  IS  WHITE  NOISE,  THE  LATTICE  IS  EQUIVALENT  TO  THE  INVERSE 
OF  THE  SYSTEM. 

***  VARIABLE  DECLARATIONS**** 

DSEED  =  SEED  USED  BY  THE  IMSL  WHITE  GAUSSIAN  NOISE  FUNCTION  GGNQF 
GRAFP  =  MAX  X-AXIS  VALUE  FOR  X  AND  Y  GRAPHS 
GRAFN  =  MAX  X-AXIS  VALUE  FOR  MSE  GRAPH 

H  =  ARRAY  OF  AUTOREGRESSIVE  PARAMETERS  THAT  DEFINE  THE  SYSTEM 
IY  =  SEED  USED  BY  THE  UNIFORM  WHITE  NOISE  FUNCTION  URAND 
MN  =  N  *  N 

MNP1  =  MN  +  1:  NUMBER  OF  ROWS  IN  THE  NONLINEAR  LATTICE 
MSE  =  MEAN  SQDaRE  ERROR  BETWEEN  THE  INPUT  NOISE  SIGNAL  AND  THE 
FORWARD  ERROR  SIGNAL  FROM  THE  TOP  ROW  OF  THE  LATTICE 
MSESUM  =  RUNNING  TOTAL  USED  IN  CALCULATION  MSE 
MSEMAX  =  MAXIMAUM  VALUE  OF  MSE  FOR  USE  IN  PLOTTING 
MSEPLT  =  REAL*4  VALUES  OF  MSE  USED  IN  DISSPLA  PLOTTING 
NORM  =  ARRAY  OF  FACTORS  THAT  NORMALIZE  THE  LATTICE  INPUTS 
N  =  DIMENSION  OF  THE  SQUARE  Y  TENSOR  MATRIX 
NUMPTS  =  NUMBER  OF  POINTS  USED  IN  CALCULATING  THE  RHO  MATRICES 
NINDEX  =  TIME  INDEX  ARRAY  USED  IN  DISSPLA  PLOTS 
PLTPTS  =  NUMBER  OF  DATA  POINTS  USED  IN  DISSPLA  PLOTS 
RANGE  =  +7-  RANGE  OF  UNIFORM  WHITE  NOISE 
RHO  =  ARRAY  OF  REFLECTION  COEFFICIENTS 
RHOSUM  =  ARRAY  USED  IN  CALCULATING  THE  AVERAGE  RHO  MATRIX 
SCALE  =  RECIPROCAL  OF  THE  NORM  OF  THE  SYSTEM  OUTPUT  FOR  WHEN 
THE  SYSTEM  IS  EXCITED  BY  WHITE  NOISE 
STDDEV  =  STANDARD  DEVIATION  OF  GUASSIAN  WHITE  NOISE 
TRIAL  =  NUMBER  OF  NOISE  REALIZATIONS  USED  IN  CALCULATING  AVG  RHO 
X  =  INPUT  INTO  SYSTEM  H(Z) 

XHAT  =  LATTICE  OUTPUT  THAT  APPROXIMATES  X 

XHPLOT  =  ARRAY  OF  REAL*4  VALUES  OF  XHAT  USED  IN  DISSPLA  PLOTS 
XMAX.XMIN  =  RANGE  OF  X  AND  XHPLOT  VALUES  USED  IN  DISSPLA  PLOTS 

y  =  Output  of  h(Z) 

YPLOT  =  ARRAY  OF  REAL*4  VALUES  OF  Y  USED  IN  DISSPLA  PLOTS 
YMAX.YMIN  =  RANGE  OF  YPLOT  VALUES  USED  IN  DISSPLA  PLOTS 
YSUM  =  RUNNING  TOTAL  USED  IN  CALCULATING  THE  SYSTEM  OUTPUT 


***VARIABLE  DECLARATIONS**** 

INTEGER  TRIAL. NUMPTS, PLTPTS, IYJM1.JM1 
REAL*4  X (5000 ) ,MSEMA$, RANGE jSTODEV’YMA 
REAL*4  XHPLOT(OOOO) , NiNDEXf $000) .MOEPL 
PEAL*8  Y( 5000) , R( 26 ,26) , RH0SUM( 20,26), 


MSEMA# . RANGE ’ STODEV ’ YMAX , YMI N . XMAX , XM I N . GRAFN 


*.  J1.  .'_H  1 


oooooooooo  cn  oo  o<~>  ro>-*  ooooo  c~i  nsien  ooo 


SET  THE  NUMBER  OF  WHITE  NOISE  REALIZATIONS  USED  TO  CALCULATE 
THE  AVERAGE  RHO  MATRIX 
TRIAL  =  25 

DEFINE  MODEL  PARAMETERS 
N  =  2 

MN  =  N  *  N 
MNP1  =  MN  +  1 
DEFINE  H(Z )  COEFFICIENTS 

THE  AR  PARAMETERS  H(2,l)=0.6  AND  H(l,2)=. 08  CORRESPOND  TO 
GAMMA(2)=-. 55  AND  GAMMA(3)=-.  08 
DO  7  1=1,  N 
DO  6  J=1.N 

H(i,j)  =  0. 

CONTINUE 
CONTINUE 
H(l.l)  =  0. 5 
HfM)  =  0.6 
H(  j  J)  =  0.  2 
HC  1.2)  =  0.08 
H(2,2)  =  0.2 
H(2  3)  =  0.  5 
H(  3  ’  1)  =  0.2 

INITIALIZE  MSESUM,  MSEMAX  AND  RHOSUM  MATRIX 
DO  2  1=1 , MNP 1 
DO  1  J=1AMNP1 

RHdSUM(I,J)  =  0. 

CONTINUE 
CONTINUE 
MSESUM  =  0. 

MSEMAX  =  0. 

XMAX  =  0. 

YMAX  =  0. 

DEFINE  SEED  VALUESa  SATURATION  LIMIT.  RANGE  OF  UNIFORM  NOISE. 
STD  DEV  OF  GAUSS  NdlSE,  AND  NUMBER  Of  POINTS  IN  TIME  SEQUENCES 
IY  =  1354 

DSEED  =  1243073. 5D0 

SAT  =  0.  7 

RANGE  =  1.73205- 

STDDEV  =  1.  0 

NUMPTS  =  5000 

PLTPTS  =  5000 

GRAFN  =  FLOAT? NUMPTS  +  1) 

GRAFP  =  FLOAT? PLTPTS  +  1) 

WRITE  HEADER  FOR  UNIFORM  NOISE  INPUT 
WRITE(8.4)  RANGE 

FORMATCTi /INPUT  WHITE  UNIFORM  NOISE  HAS  ZERO  MEAN  AND  RANGE  OF 
*+/-  1  F6. 3) 

8.5) 

TT  •  D 


WRITE? 8.5;  TRIAL, NUMPTS. IY  ,  , 

FORMATfTl/RHO  IS  AVERAGED  OVER  ',13/  TRIALS  OF  ',15/  POINTS. 

*  INITIAL  SEED-  ,16) 

WRITE  HEADER  FOR  GUASSIAN  NOISE  INPUT 
WRITEC 8 . 4)  STDDEV 

FORMAT? tl,'  INPUT  WHITE  GAUSS  NOISE  HAS  ZERO  MEAN  AND  STD  DEV  OF1 
*F6. 3) 


WRITEC8 , 5)  TRIAL.NUMPTS, DSEED 

FORMAT? tl,  RHO  1$  AVERAGED  OVER  ',13/  TRIALS  OF  ',15,'  POINTS. 
*  INITIAL  SEED=f,F10. 1) 

PRINT  H  MATRIX 
WRITE? 8. 8) 

FORMATCTf/Y(K)=X(K)-(TENSOR  PRODUCT  H*Y)  WHERE  H(  I ,  J)  =  ') 

DO  10  1=1, N 

WRITE(8 ,9)  (H(I,J),J=1,N) 
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FORMAT (T10,5(F6. 3,4X)) 

CONTINUE 

** r************************************************************** 

RUN  TRIALS  FOR  DIFFERENT  SEED  VALUES  TO  GET  AVERAGE  RHC  MATRIX 
****************************** 

DO  50  L=l, TRIAL 

DSEED  =  DSEED/DFLOAT( L) 

IY  =  IY/L 

SELECT  EITHER  A  GUASSIAN  OR  UNIFORM  INPUT  NOISE  SEQUENCE 
DO  11  K=1 .NUMPTS 

THE  INPUT  NOISE  IS  UNIFORM  ON  (-RANGE , RANGE)  WITH 
MEAN  =  0  AND  VARIANCE  =((2*RANGE)**2)/'12) 

X(K)  =  2.0  *  (URAND(Iy|  -  .5)  *  RANGE 
THE  INPUT  NOISE  IS  GAUSSIAN,  MEAN=0  AND  VARIANCE=STDDEV**2 
X(K)  =  GGNQF(DSEED)  *  STDDEtf 
CONTINUE 


mxm :  m 


)  +  H(2, 1)*Y( 1)  +  H(3,1)*Y(1)*Y(1)) 


YSUM  =  0. 

DO  14  1=1. N 

DO  i3  J=1 ,N 
IM1  =  1-1 
JM1  =  J-l 

YSUM=H(I ,J)*C00RD(Y(K-1) , I M 1 ) * COORD( Y ( K-2 ) ,JM1)+YSUM 
CONTINUE 
CONTINUE 

Y(K)  =  D6LE(X(K))  -  YSUM 
CONTINUE 

***CALL  SUBROUTINES**** 

DETERMINE  AUTOCORRELATION  MATRIX  FOR  Y  SEQUENCE 
CALL  NLCLAT(Y  NUMPTS .R, N) 

DETERMINE  REFLECTION  EaUtORS  FROM  AUTOCORRELATION  MATRIX 
CALL  SCHURf RHO. R, ALPHA, BETA. MNP 1 ) 

ADD  TOGETHER  THE  AHO  MAfRICE^  FROM  EACH  TRIAL 
DO  45  1=1 ,MNP1 
DO  40  J=1aMNP1 

RHO§UM(I,J)  =  RHO(I.J)  +  RHOSUM(I.J) 

CONTINUE 

CONTINUE 

CONTINUE 

CALCULATE  AVERAGE  RHO  MATRIX  AND  TRUNCATE  TO  TWO  DECIMALS 
DO  54  1=1 ,MNP1 
DO  53  J=1 ,MNP1 

RHO(I.J)  =  RHOSUM(I,J)/DFLOAT(TRIAL) 

RHO(I,J)  =  OINT(RHO( I , J)*100.  )/100. 

CONTINUE 

CONTINUE 

WRITE! 8, 55) 

FORMAT(Ti,  RHO( I .J)  =') 

DO  60  I  =  l.MNPl 

WRITE (8 , 58)( RHO( I , J) , J=1 ,MNP1 ) 

FORMAT(T1,10F6.  2) 

CONTINUE 


***  NORMALIZE  THE  LATTICE  INPUT  SIGNALS  AND  PASS  THE  *** 

NOISE  GENERATED  DATA  THROUGH  THE  LATTICE  FILTER. 

CALL  NORMS(Y,NANUMPTS,NORM) 

SCALE  =  1. 0/NOftM(l) 

CALL  NLLAT(Y, RHO ,N, NUMPTS, NORM, XHAT) 

***EXAMINE  THE  WHITENING  EFFECT  OF  THE  LATTICE  FILTER  BY  ***** 
CALCULATING  THE  MEAN-SOARE  ERROR  BETWEEN  THE  INPUT  WHITE 
NOISE  AND  THE  FORWARD  ERROR  SIGNAL  OUTPUT  FROM  THE  TOP 
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P. v 


K. 


63 

C 

C 

c 

c 

c 


85 

C 


88 

89 


90 

C 

C 

C 

C 

C 

C 

C 


98 


99 

C 

C 

C 

C 

C 


ROW  OF  THE  LATTICE. 

DO  63  K=1 .NUMPTS 

MSESUPf  =  ^DBJ_E^X^K^-XHAT^K^*rDBLE(X(K))-XHAT(K))  +  MSESUM 

mse[>l)(k§  =Rsngumse(k))°  T  ,  v 

IF(_MSEPLT(K).  GT.  MSEMAX;  MSEMAX  =  MSEPLT(K) 

CONTINUE 


NOW  THAT  THE  INVERSE  FILTER  HAS  BEEN  EVALUATED,  DO  THE 
DECONVOLUTION  PROBLEM.  FIRST.  DEFINE  THE  "UNKNOWN"  INPUT 
SEQUENCE  X(K)  AND  THEN  GENERATE  THE  SYSTEM  OUTPUT  Y(K). 


DO  85  K=1 .PLTPTS 

IF(K.  LE.250)  X!K)  =  FLOAT ( K)/250. 
IF(K.  LE.  750.  AND.  K.  GT.  250 )  w/" ' 

I  F(  K.  LE.  1250.  AND.  K.  GT.  750) 

I  F(K.  LE.  1750.  AND.  K.  GT.  1250 
IF!K.  LE.  2000.  AND.  K.  GT.  1750 


IF 

IF 

IF 

IF 

IF 

IF 

CONTINU 


K.  LE.  2400.  AND.  K.  GT.  2000 
K.  LE.  2700.  AND.  K.  GT.  2400 
K.  LE.  2900.  AND.  K.  GT.  2700 
K.  LE.  3000.  AND.  K.  GT.  2900 
K.  LE.  3500.  AND.  K.  GT.  3000 
K.  LE.  5000.  AND.  K.  GT.  3500 


X(K)  =  2.0  -  FLOAT!K)/250. 

X!  K)  =  FLOAT ( K 1/ 250.  -  4.0 
X ( K )  =  1.  0  -  FLOAT! K- 1250) /2 50. 
FLOAT ( K-1750)/250.  -  1.  0 

-0.  5 
0.  5 
-0.  5 

0.  7 * S I N ( 0 .  0126* FLOAT ( K— 3001 ) ) 


(9= 


0. 9*SIN(0. 0021* FLOAT ( K— 350 1 ’ 


}j) :  m 


)  +  H( 2 , 1 )*Y( 1 )  +  H( 3 , 1 )*Y( 1 )*Y( 1 ) ) 


Y( 1 )  =  DBLE(X(i; 

Y!2  j  =  DBLE( X{ 2 
DO  90  K=3 ,  PLTPTI 
YSUM  =  0. 

DO  89  1=1 .N 

DO  98  J=1,N 
IM1  =  1-1 
JM1  =  J-l 

YSUM=H( I ,  J)*COORD( Y( K- 1 ) , IM1 )*COORD(Y( K-2) ,  JM1 )+YSUM 
CONTINUE 
CONTINUE 

Y(K)  =  DBLE(X(K))  -  YSUM 
CONTINUE 


NOW  PASS  THE  SYSTEM  OUTPUT  DATA  THROUGH  THE  LATTICE  FILTER 
EMBEDDED  WITH  THE  PREVIOUSLY  CALCULATED  REFLECTION  FACTORS 
RHO^I  J).  THE  FORWARD  ERROR  SIGNAL  OUT  OF  THE  TOP  ROW  OF  THE 


LATl 


IS  XHAT^K^  AND. SHOULD. APPROXIMATE  THE  DESIRED  SIGNAL 


AFTER  IT 


^NORMALIZED  AND  SCALED. 


CALL  NORMS! Y , N .PLTPTS .NORM) 

CALL  NLLATfY  RHQ.N .PLTPTS , NORM.XHAT) 
DO  98  K=1 , PLtPTS 


XHAT!^)  =  SCALE  *  (XHAT(K)*NORM(  1 ) ) 
XHPLOt(K)  =  SNGLfXHAT ( K) ) 

YPLOT(K)  =  SNGL(YfK)) 

IF(ABS(X(K)). GT. XMAX)  XMAX  =  ABS(X(K 
IF(ABS(XHPLOT(Kj). GT. XMAX)  XMAX  =  AB 


I  F(_ABS(  YPLOT(K) 
CONTINUE 
XMIN  =  -1.0*  XMAX 
YMIN  =  -1.0  *  YMAX 
DO  99  K=1  5000 

NINDE)C(K)  =  FLOAT(K-l) 
CONTINUE 


ll 


XHPLOT 


GT.  YMAX)  YMAX  =  ABS(YPLOT(K 


«» 


DISSPLA  PLOTTING  ROUTINES 


PLOT  SYSTEM  INPUT  AND  LATTICE  OUTPUT 
CALL  TEK618 
CALL  COMPRS 
CALL  RESET! 'ALL') 

CALL  PAGE(8.  0,6.  5) 
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LI  ,1  iLJWUWIH.  — J 


c 

c 


CALL  XINTAX 

CALL  AREA2D(6.  75,5.  0) 

CALL  XNAMEf'TIME  SAMPLES* ' , 100) 

CALL  YNAMEf ’X(KJ,XHAT(K)$' ,100) 

CALL  YNAMEpX(K)v  ,  100) 

CALL  CROSS 

CALL  GRAFfO. ,500, ,GRAFP ,XMIN , 1  SCALE' ,XMAX) 

CALL  LINESP  (2.  OV 

CALL  LINESfXfK}*'  IPAK.l) 

CALL  LINESf 1 XHATC K J$ 1 ,lPAK,2) 

CALL  LEGLIN 
CALL  DOT 

CALL  CURVEfNINDEX ,X , PLTPTS ,0) 

CALL  RESET  'DOTV 

CALL  CURVEfNINDEX ,XHPLOT. PLTPTS ,0) 

CALL  LEGEND  (IPAK  2,5.  0,(5.  5) 

CALL  ENDPL(O) 


PLOT  Y 


CALL  AREA2D(6.  75,5.  0) 

CALL  XNAMEpTIME  SAMPLESS' ,100) 

CALL  YNAMEf 'Y(K)$r, 100) 

CALL  HEADINC'YCK)  =  X(K)  -  0. 2*Y(K-1)**2$' .100,1.  5.1) 

CALL  HEADINf 'Y(K)  =  XfKW  0.  6*Y[  K-1>0.  08VY?  KL2) W '  100 . 1,  5 , 1) 
CALL  HEADINf  ^Y(K)=X(K)-(0.  1*Y(K-I)+  0. 5*Y( K-I)*Y( K-2 1 , 


*100,1. 5.1 , 

caLl  cPoss 


CALL  GRAFfO. ,500.  ,GRAFP ,YMIN .  SCALE  , YMAX ) 
CALL  CURVEfNfNDEX  YPLOT  PLTptS.O) 

CALL  ENDPL(O) 


PLOT  MSE 

CALL  AREA2D(6.  75,5.  0) 

CALL  XNAMEPTIME  SAMPLESS'.IOO) 

CALL  YNAMEf 1 MSE(K)$ 1  , 100) 

CALL  GRAFfO. ,500. ,GrAFN,0. , 'SCALE' ,MSEMAX) 
CALL  CUR VhfN INDEX ,MSEPLt,NUMPTS ,0) 

CALL  ENDPLfO) 

CALL  DONEPl 


STOP 

END 

£*********************************************************************** 

c 

C  SUBROUTINE  SCHUR 
C 

C  CALCULATES  THE  REFLECTION  FACTORS  FROM  THE  CORRELATION  MATRIX 
C 

C  WRITTEN  29  APRIL  1985  BY  P.J.  LENK  (REF  12: P.  174) 

Q*********************************************************************** 

SUBROUTINE  SCHURfRHO.R .ALPHA, BETA, N) 

REAL*8  RHOf  26 , 26 )  ,R(26,26)  ,Af.PHA(26 ,26) , BETAf  26 , 26 )  ,RNORM,T 


INITIALIZE  THE  ALPHA  AND  BETA  ARRAYS 


DO  10  I  =  1,N 
DO  5  J  = l.N 


5 

C 

C7 

10 

C 

C 

c 


®S^L-4SAfKSr(R(I’I)) 

RHOf  1,0)  =  0.0 


CONTINUE 
CONTINUE 


BEGIN  CALCULATING  THE  REFLECTION  FACTORS 
DO  50  J  =  2 ,N 
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Lv 

A  A  » -*  A  Ivl'n  .v  A  A.  A  .'-A-a' 


nnnn  ooo  ooo  c->c->o  ooo  ooooooooo 


30 

40 

C 

C42 

C 

50 

C 

C 


NJ1  =  N-  J  +  1 
DO  40  I  =  1 , N J 1 
JI1  =  J  +  I  -  1 
IP1  =  I  +  1 


RH0( I , J 1 1 )  =  ALPHA( I , JI 1)/BETA(IP1 , JI1 ) 
RNORM  =  DSQRT( 1.  0  -  ftH0(I ,JI1)*RH0( I ,Jll)) 


DO  30  K  =  1,N 
T  =  ALPtfAfl.K 


ALPHA! I . K)  =  f ALPHAfI ,K)-RH0( I ,JI1 

n  f--r  i  /  V  A  \  '  /  (Ar-r  1  /  T  p  |  J  m  in  >  ▼  '  i  -r  i 


BETA( I , d)  =  (BETA(I 
CONTINUE 
CONTINUE 


RHO(i;JIl)*T)/R 


BETA^ 


CONTINUE 


WR IT E(8 ,42)J,( (ALPHA! I,K).K=1,N),I=1, 

FORMAT  (Y2X, l3,4!2X. El2. 5)1 

WRITE( 8 , 42 ) J , ( ( BETA( I , K) , K=1 , N) , 1=1 , N 


=1,N) 
N) 


IP1 ,K))/RN0RM 
ORM 


RETURN 

END 


*********************************************************************** 


FUNCTION  URAND 

TAKEN  FROM  "COMPUTER  METHODS  FOR  MATHEMATICAL  COMPUTATIONS"  BY 
G.  E.  FORSYTHE,  M.  A.  MALCOLM,  AND  C.  B.  M0LER 

*********************************************************************** 


C 

C 

C 

c 

c 


REAL  FUNCTION  URAND(IY) 

INTEGER  I A  I C , I  TWO , M2 ,M,MIC 
DOUBLE  PRECISION  hALFM 
REAL  S 

DOUBLE  PRECISION  DATAN.DSQRT 
DATA  M2/0/.ITWO/2/ 

I F(M2. NE.  0)  GO  TO  20 

IF  FIRST  ENTRY,  COMPUTE  MACHINE  INTEGER  WORD  LENGTH 


M=1 
10  M2=M 

M=ITW0*M2 

IF(M.  GT.M2)  GO  TO  10 
HALFM=M2 

COMPUT  MULTIPLIER  AND  INCREMENT  FOR  LINEAR  CONGRUENTIAL  METHOD 

IA=8*IDINT(HALF,1*DATAN(  1.  D0)/8.  D0)  +  5 
IC=2*IDINT( HALFM*(. 5D0-DSQRT(3.  D0)/6.  DO) )  +  l 
MIC=(M2-IC)+M2 

S  IS  THE  SCALE  FACTOR  FOR  CONVERTING  TO  FLOATING  POINT 
S=. 5/HALFM 

COMPUTE  NEXT  RANDOM  NUMBER 
20  IY=IY*IA 

THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHICH  DO  NOT  ALLOW 
INTEGER  OVERFLOW  ON  ADDITION 

IF(IY.GT.MIC)  IY=( IY-M2)-M2 

IY=IY+IC 


THE  FOLLOWING  IS  FOR  COMPUTERS  FOR  WHICH  THE  WORD  LENGTH 
FOR  ADDITION  IS  GREATER  THAN  FOR  MULTIPLICATION 
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»  mmiuiiii 


j  w  y  M  V  » 


IF( IY/2. GT.  M2)IY=( IY-M2)-M2 

HE  FOLLOWING  STATEMENT  IS 
■FFECTS  THE  SIGN  BIT 

I  F(  I Y.  LT. 0)  IY=(  IY+M2)+M2 
"3AN0=FL0AT(:' 


THE  FOLLOWING  STATEMENT  IS  FOR  COMPUTERS  WHERE  INTEGER  OVERFLOW 
AFFECTS  THE  SIGN  BIT 


IY)*S 


UR 

RETURN 
END 

p********************* A ************************************************* 

C  * 

C  SUBROUTINE  NLCLAT 

C  * 

C  THIS  SUBROUTINE  PRODUCES  A  CORRELATION  MATRIX  FROM  NONLINEAR  * 

C  TIME  SEQUENCE  IN  AN  ORDER  WHICH  IS  COMPATIBLE  WITH  SUBROUTINE  * 

C  SCHUR.  * 

C  * 

C  WRITTEN  7  MAY  1985  BY  P.J.  LENK  (REF  12: P.  181) 

C 

£************************************************************★********** 


10 

20 

C 

C 

C 


30 


40 

50 

C 

C 

C 


C 

C12 

60 

70 

80 

C 


SUBROUTINE  NLCLAT  (Y.IYS.R.N) 
REAL*8  Y(5000),R(26,26),$UM,V 

DEFINE  CONSTANTS 


EC( 26) 


MN  =  N*N 
MNP1  =  MN  +  1 
IYSM2  =  IYS  -  2 
FIYSM2  =  FL0AT( IYSM2) 


INITIALIZE  R  MATRIX  TO  ZERO 


DO  20  I  =  l.MNPl 
DO  10  J  =  l.MNPl 

r\  /  t  i  \  1  r\  r\ 


R(  I ,  J)  =  0.0 
CONTINUE 


CONTINUE 
BEGIN  OUTER  LOOP 


DO  80  I  ~  3, IYS 
IR  =  1 

VEC(IR)  =  Y( I) 

DO  50  MP1  =  1 ,  N 
MO  =  MP1  -  1 
LLIM  =  2*MP1  -  1 
DO  40  L  =  1 , LLIM 
LO  =  L  -  1 
II  =  MO 
J1  =  LO/2 

IF  (MOD( LO ,2). EQ. 0)  GO  TO  30 
II  =  J1 
J1  =  MO 
IR  =  IR  +  1 

VEC(IR)  =  COORD(Y( 1-1) , 1 1 ) *COORD( Y( 1-2) ,J1) 
CONTINUE 
CONTINUE 


CALCULATE  THE  CORRELATIONS 


DO  70  J  =  l.MNPl 
DO  60  K  =  J.MNPl 


CONTINUE 

CONTINUE 

CONTINUE 


R(J , K)  =  R(  J , K)  +  VEC( J)*VEC( K) 
WRlTE[6p2)VEC(J),VEC(K),R(J,K) 
FORMAT (3(2X,E12.  5)) 
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|y>~y'VLV„\-  -'.V.v.v.v.v/.-  «VT J"*- V.'-  L-". 

s.'s  v_  vO  s  *.-  V-V-V'-  V- vLV- <*•- V-V- «’.VuVlVv  ^-V'-V'.V-  i.“-V iv.  >*.  tf/M 


c 

c 


DIVIDE  BY  THE  NUMBER  OF  DATA  ELEMENTS  CONSIDERED 


DO  100  J  =  1 , MNP  1 
DO  90  K  =  J , MNP 1 

R( J , K)  =  R( J , K)/FIYSM2 
90  CONTINUE 

100  CONTINUE 
C 

C  FILL  IN  THE  SYMMETRIC  HALF  OF  CORRELATION  MATRIX 
C 

DO  120  I  =  2 ,MNP1 
IM1  =  I  -  1 
DO  110  J  =  1 , IM1 
R( I , J)  =  R(J,I) 

110  CONTINUE 

120  CONTINUE 
C 
C 

RETURN 

END 

C 

C 

£*********************************************************************** 

C 

C  FUNCTION  COORD 

C 

C  GENERATES  OUTPUT  OF  RANDOM  FUNCTION 

C  CREATED  23  AUG  84  (REF  12: P.  187) 

£***************** ******************************************************* 
DOUBLE  PRECISION  FUNCTION  COORD(X.I) 

C  USE  SIMPLE  POWER  SERIES  TYPE  POLYNOMIALS 

C 

Y  =  1.0 

IF  (I. Eg. 0)  GO  TO  jO 
y  =  X**I 

30  COORD  =  DBLE(Y) 

RETURN 

END 

C 

c 

£***** ******************************************** ******************** 

C 

C  SUBROUTINE  NLLAT 

C 

C  THIS  SUBROUTINE  IMPLEMENTS  THE  NONLINEAR  LATTICE  FILTER  USING 
C  PREVIOUSLY  CALCULATED  REFLECTION  FACTORS. 

C 

C  WRITTEN  30  APRIL  86 

C 

Q* *************************************************************** ******* 

SUBROUTINE  NLLATf Y , RHO  ,N .NUMPTS ,NORM,XHAT) 

C  ***VARIABLE  DEFINITIONS^*1** 

C 

C  INFWD(I)  =  FORWARD  ERROR  INPUT  INTO  THE  1 1 TH  ROW  OF  THE  LATTICE 
C  INBKDfl  =  BACKWARD  "  r'  11  ''  "  11 

C  OUTFWD? I )=  FORWARD  ERROR  OUTPUT  FROM  THE  1 1 TH  ROW  OF  THE  LATTICE 
C  OUTBKDC  n=  BACKWARD  "  "  "  ‘‘  "  "  •'  h  " 

C  Y  =  INPUT  DATA  VECTOR 

C  RHO  =  REFLECTION  FACTOR  MATRIX 

C  NORM  =  VECTOR  OF  NORMS  OF  THE  LATTICE  INPUT  TERMS 

C  XHAT  =  OUTPUT  DATA  VECTOR;  IT  IS  THE  FORWARD  ERROR  SIGNAL  FROM 

C  THE  LAST  STAGE  OF  THE  FIRST  ROW 

C  NUMPTS  =  NUMBER  OF  POINTS  IN  THE  INPUT/OUTPUT  SEQUENCES 

C  N  =  DIMENSION  OF  SQUARE  Y  DATA  MATRIX 

C  MNP1  =  DIMENSION  OF  THE  RHO,  NORM,  AND  INPUT/OUTPUT  ARRAYS 

C 

C  ***VARIABLE  DECLARATIONS**** 

INTEGER  N,MN, NUMPTS, LAST, MNP1 


y y-i AAstfOV  A.- .-v 


ffVl  fftPTl'VT'TW 


^1'  v 


WJW  ‘j-w  7  ■  T 


c 

c 


10 


11 

12 

13 

C 

C 


14 

C 

C 

15 


C 

C 

20 


REAL*8  Y(5Q00) .RHO(26 ,26) ,NORM(26),XHAT(5000) , INFWD(26) 

REAL*8  I NBKD( 26 )  ,OUTFWD(26)  ,0UTBKD(26)  ,RN0RM 

MN  =  N  *  N 
MNP1  =  MN  +  1 

INITIALIZE  THE  INPUT  AND  OUTPUT  ARRAYS  TO  ZERO 
DO  5  1=1.26 

INFWCffl)  =  0. 

INBKDfl)  =  0. 

OUTFWDf  I )  =  0. 

OUTBKD(I)  =  0. 

CONTINUE 

DETERMINE  THE  LATTICE  INPUTS  FOR  EACH  TIME  K. 

DO  60  K=1 .NUMPTS 

CLEAR  INPUTS  FOR  THIS  TIME  ITERATION 
DO  10  1  =  1 , MNP  j 

INFWD(I)  =  0. 

INBKD(I)  =  0. 

CONTINUE 

I  F(  K.  EO.  1)  GO  TO  15 
IF(K-  EO.  2)  GO  TO  20 
SET  LATTICE  INPUTS  FOR  CASE  WHEN  K>=3 
IR  =  1 

INFWD(IR)  =  Y( K)/NQRM( IR) 

DO  13  MPl  =  i:n 

MO  =  MPi  -  1 
LLIM  =  2*MP1  -  1 
DO  12  L  =  1, LLIM 
LO  =  L  -  1 
II  =  MO 
J1  =  LO/2 

IF  (MOD( L0,2). EQ. 0)  GO  TO  11 
II  =  J1 
J1  =  MO 
IR  =  IR  +  1 

INFWD(IR)=COORD(Y(K-l),Il)*COORD(Y(K-2),Jl)/NORM(IR) 


CONTINUE 
CONTINUE 
INFWD(1 


;d 


c 

INFWD 

3 

K-l 

/NOkM(3 

c 

INFWD 

4 

K-2 

/N0RMC4 

c 

INFWD 

5 

K-l 

*Y(K-2) 

c 

INFWD 

6 

K-l 

*Y( K-l ) 

c 

INFWD 

7 

K-2 

|*Y( K-2 1 

c 

INFWD 

8 

,K-1 

)*Y( K-lj 

c 

INFWD 

9 

=  Y 

K-l 

|*YC K-2) 

c 

INFWD 

DO 

*  id 
[4 

))  = 
J=l> 

tr 

/NORM 

/NORM 

/NORM 

*Y(K-_ 

*Y(K-2 

)  *  Y  ( K 


/N0RMT8 
2)/NORM(9 


Y(  K-2)/NORM( 10) 


INBKD(J)  =  INFWD(«J+1) 

CONTINUE 
GO  TO  25 

SET  INPUTS  FOR  K=1  CASE 
INFWD(l)  =  Y( 1  )/NORM( 1 ) 
INFWDT2)  =  l./NORM(2) 
INBKD(l)  =  INFWD(2) 

GO  TO  25 


SE 

IN 


:WD(  1 
NFWD 
NFWD 
NFWD 
NFWD 
NFWD 
NBKD 
NBKD 


rs  FO  _  .  . 

=  Y( 2)/NORM( 1 

2)  =  1./NORMC2 

3)  =  Y( 1 )/NORM 
6)  =  Y  *Y1J 


!): 


i*Y( 

1) '=  INFWDT2 

2)  =  INFWD(3 


Y^1)/N0RM( 18) 


/R9RMQ1J 
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l 


INBK0(5)  =  INFWD(6) 

INBKDflO)  =  INFWD(ll) 

INBKDC 17)  =  INFWD( 18) 

IMPLEMENT  THE  J'TH  UPPER  DIAGONAL  OF  THE  LATTICE.  MOVE  ALONG  THE 
DIAGONAL  BY  STARTING  AT  ROW  L=1  AND  MOVING  DOWN  TO  ROW  L=LAST. 

DO  50  J=1.MN 

LXST  =  MNP1  -  J 
IFfJ.EQ.  n  GO  TO  40 

IF_NOT  ON  THE  FIRST  UPPER  DIAGONAL  (I.E.  J  NOT  =  1)  THEN  UPDA 


UPPER  DIAGONAL  [I.E.  J  NOT  =  1)  THEN  UPDATE 
OUTPUTS  OF  THE  PREVIOUS  DIAGONAL 


THE  INPUTS  FROM  THE  OUTPUTS  OF  THE  PREVIOUS  DIAGON, 
DO  35  L=1 , LAST 

INFWD(L)  =  OUTFWDf  L) 

INBKDC L)  =  OUTBKD  L+l) 

CONTINUE 

DO  THE  LATTICE  CALCULATIONS 
DO  45  L=1 , LAST 


(  L?S=R[|NFWDf  L?H-(  RflO?  L?J  WlftBKD?  f }  V  RNORM 
(L)  =  (INBKDC L)  -  RHO(L!vJ+L)*INFWD(L))/RNORM 


60  CONTINUE 
RETURN 
END 


CONTINUE0™5  =  (INBKD(l5  "  RHOCL,J+L)*INFWD(L) ;/ 
CONTINUE 

REMOVE^  NORMALIZATION^  FACTOR  FROM  THE  OUTPUT  OF  THE  TOP  ROW 


C********************************************************************** 

C 

C  SUBROUTINE  NORMS 

l  THE  N0RHS  0F  THE INPUTS 

c  WRITTEN  30  APRIL  86 

C 

£*********************************************************************** 
SUBROUTINE  NORMS(Y,N,NUMPTS,NORM) 


1 LLIM' 1 L0  ■ 'M0  ’ 'HP1  ■ 11  ■ J1  ■ J - L 

MN  =  N  *  N 
MNP1  =  MN  +  1 

C  INITIALIZE  VECTORS  TO  ZERO 

DO  10  K  =  1 , MNP  1 

NOR©)~=°b. 

10  CONTINUE 
C 

C  SINCE  TIME  INDEX  STARTS  AT  1=3.  INITIALIZE  NORMS  OF  INPUTS  WHICH 

C  ARNEORPM  1  f=  °Y  1  (I*'  1^°+T2^  5?8)™ES  1=1  AND  I=2‘ 

NORM  2)  =  1.0  +1.0 


NORMf 3 )  =  Y[l)  *  Y( 1 ) 

NORMf 6)  =  NORMC 3)  *  N0RM(3) 
NORMf 11)  =  NORM(b)  *  NORMf 3) 
N0RM(18)  =  NORMf 6)  *  N0RM(6) 


WH0SE  components  match  the 

DO  80  I  =  3.NUMPTS 
IR  =  1 

VEC(IR)  =  Y( I ) 


V  /v 


30 

40 

50 

C 

C 


60 

80 

C 

83 


86 

85 


DO  50  MP1  =  l.N 

MO  =  MPi  -  1 
LLIM  =  2*MP1  -  1 
DO  40  L  =  1, LLIM 
LO  =  L  -  1 
II  =  MO 
J1  =  LO/2 

IF  (MOD( LO , 2). EQ. 0)  GO  TO  30 
II  =  J1 
J1  =  MO 
IR  =  IR  +  1 

VEC(IR)  =  COORD(Y( 1-1) ,  I1)*C00RD(Y( 1-2) , J 1 ) 

CONTINUE 

CONTINUE 

CALCULATE  THE  NORMS 
DO  60  K=1 ,MNP1 

IF(K.  E6.  2)  GO  TO  60 

NORM(K)  =  NORM(K)  +  VEC(K)*VEC(K) 

CONTINUE 

CONTINUE 

WRITE(8.83) 

FORMAT CTj'R' .T10, 'NORM(K)' ) 

DO  85  K=l,MNf)l 

N0RM(1^]  =  DSQRT(NORM(K)/DFLOAT(NUMPTS)) 

WRITE(8,86)  K.NORMCK) 

F0RMAT(Yi,I3,T8,E12.  5) 

CONTINUE 

RETURN 

END 
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